1{ 2 "cells": [ 3 { 4 "cell_type": "markdown", 5 "id": "eac3271e", 6 "metadata": {}, 7 "source": [ 8 "# The wavy ramp\n", 9 "\n", 10 "In this example, we will integrate the dynamics of a point mass bouncing down a non-linear ramp. The point mass behaves like a rubber ball, meaning that at every bounce it will dissipate part of its kinetic energy due to the [inelastic](https://en.wikipedia.org/wiki/Inelastic_collision) character of the collisions.\n", 11 "\n", 12 "Let us begin, as usual, with the definition of the dynamics:" 13 ] 14 }, 15 { 16 "cell_type": "code", 17 "execution_count": 1, 18 "id": "553b9d9f", 19 "metadata": {}, 20 "outputs": [], 21 "source": [ 22 "import heyoka as hy\n", 23 "\n", 24 "x, y, vx, vy = hy.make_vars(\"x\", \"y\", \"vx\", \"vy\")\n", 25 "\n", 26 "eqns = [(x, vx),\n", 27 " (y, vy),\n", 28 " (vx, hy.expression(0.)),\n", 29 " # Downwards constant acceleration.\n", 30 " (vy, hy.expression(-9.8))]" 31 ] 32 }, 33 { 34 "cell_type": "markdown", 35 "id": "cbe82079", 36 "metadata": {}, 37 "source": [ 38 "In other words, we are setting up a 2D dynamical system in which the only acceleration is due to the (constant) gravitational field.\n", 39 "\n", 40 "The contour of the ramp is defined by the non-linear equation\n", 41 "\n", 42 "$$\n", 43 "y = 1 - x + 0.05 \\cos\\left(11 \\pi x \\right).\n", 44 "$$\n", 45 "\n", 46 "Let us take a look:" 47 ] 48 }, 49 { 50 "cell_type": "code", 51 "execution_count": 2, 52 "id": "7faf59f9", 53 "metadata": {}, 54 "outputs": [ 55 { 56 "data": { 57 "image/png": "iVBORw0KGgoAAAANSUhEUgAAAhsAAAFlCAYAAABC5yqRAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAAv6klEQVR4nO3deXRUVYLH8d/NQsKSAIGwJATDGggmbAkiASEIJAEDpBFkV2QREcUdnVbHHhWx29FuFRUHBTdURG2UQWwZEXFhE5VFWkgLSEQlLDbIDrnzB+ihaZYiVNWt5fs5J+ek6r289/OeSP3y3q1bxlorAAAAX4lwHQAAAIQ2ygYAAPApygYAAPApygYAAPApygYAAPApygYAAPCpKFcnrlmzpk1NTXV1egAA4EWff/75dmtt4qm2OSsbqampWrFihavTAwAALzLGbD7dNm6jAAAAn6JsAAAAn6JsAAAAn6JsAAAAn6JsAAAAn6JsAAAAn6JsAAAAn6JsAAAAn6JsAAAAn6JsAAAAn6JsAAAAn6JsAAAAn6JsAAAAn6JsAAAAn6JsAAAAn6JsAAAAn6JsAAAAn6JsAAAAn6JsAAAAn6JsAAAAn6JsAAAAn6JsAAAAn6JsAAAAn6JsAAAAn6JsAAAAnwq5smGtdR0BAACcIOTKxl133aWcnBz98MMPrqMAAACFYNnIycnRoUOHlJ6errfeest1HAAAwl7IlY24uDgVFhaqoKBAw4YN0/333+86EgAAYS3kysav0tLSNHz4cD3yyCO67rrrXMcBACBshWzZkKTatWvrqquu0muvvaZrrrnGdRwAAMJSSJcNSapataquvPJKvfHGG5o4caLrOAAAhJ2QLxuSFB8fr8GDB+vJJ5/UlClTXMcBACCshEXZkKQaNWpo4MCBuv322/Xee++5jgMAQNgIm7IhSfXq1VOvXr00aNAgbdmyxXUcAADCQliVDUlq0aKFmjVrpp49e+rIkSOu4wAAEPLCrmxIUrdu3bRnzx6NGTPGdRQAAELeWcuGMeY5Y8w2Y8ya02w3xpjHjDHFxphVxpg23o/pXZGRkerXr59mzZqld955x3UcAABCmidXNmZIyj/D9gJJTY5/jZH01PnH8r34+HgVFBRo5MiR2rlzp+s4AACErLOWDWvtR5LO9GrcR9IL9pglkqoZY+p6K6AvXXjhhUpKStKQIUNcRwEAIGR5Y85GsqQT39pRcvy5oJCfn68lS5Zo+vTprqMAABCSvFE2zCmes6fc0ZgxxpgVxpgVpaWlXjj1+atYsaIuu+wy3XLLLdqxY4frOAAAhBxvlI0SSSknPK4naeupdrTWPmOtzbLWZiUmJnrh1N7RuHFj1a9fX1dffbXrKAAAhBxvlI23JQ0//q6U9pL+aa39wQvH9asePXrogw8+4N0pAAB4mSdvfX1F0meS0owxJcaYkcaYscaYscd3mSfpW0nFkv5H0jifpfWhypUrq0ePHho7dqwOHDjgOg4AACEj6mw7WGsHnWW7lXSd1xI5lJGRodWrV2vChAmaOnWq6zgAAISEsFxB9HSMMSooKNBLL72k1atXu44DAEBIoGycJCEhQe3atdPIkSNdRwEAICRQNk4hJydHxcXFmjFjhusoAAAEPcrGKURHRys/P1+333679u3b5zoOAABBjbJxGmlpaUpISNCNN97oOgoAAEGNsnEGeXl5evnll7V27VrXUQAACFqUjTNgsigAAOePsnEWOTk5Wr9+vV544QXXUQAACEqUjbOIjo5WXl6ebrvtNu3fv991HAAAgg5lwwNpaWmqXr26br75ZtdRAAAIOpQNDxhjlJeXpxdeeEFff/216zgAAAQVyoaHEhISdNFFF2nEiBGuowAAEFQoG+egQ4cOKi4u1vPPP+86CgAAQYOycQ5+XVn01ltv1d69e13HAQAgKFA2zlHTpk2VmJioCRMmuI4CAEBQoGyUQ15enmbOnKmvvvrKdRQAAAIeZaMcqlWrpg4dOjBZFAAAD1A2yuniiy9WSUmJ/vznP7uOAgBAQKNslFNUVJQKCwt1zz336Pvvv3cdBwCAgEXZOA/169fXhRdeqMGDB7uOAgBAwKJsnKfc3FytWbNG06ZNcx0FAICARNk4TxUqVNBll12mW2+9VT/++KPrOAAABBzKhhc0bNhQaWlp6tevn6y1ruMAABBQKBte0r17dxUXF+uBBx5wHQUAgIBC2fCS6Oho9evXTw8++KCWLVvmOg4AAAGDsuFFtWvXVm5urvr166c9e/a4jgMAQECgbHhZVlaWEhIS1KtXL5WVlbmOAwCAc5QNLzPG6LLLLtM//vEP3XDDDa7jAADgHGXDBypUqKArrrhCL7zwgp544gnXcQAAcCrKdYBQVbVqVQ0aNEh33HGHqlWrpqFDhzrN8/PPP2vu3LlatGiR1qxZo9LSUu3bt0/GGMXHx6tevXrq2LGjrrjiCjVr1sxpVgBAaKFs+FBSUpIGDBigsWPHqlKlSvrd737n1/Pv3r1bU6dO1axZs7R69WrVqlVL9erVU926ddWsWTPFxsbKWqu9e/dqx44dmj17th566CE1aNBA48eP1zXXXKOICC5+AQDOj3G1CFVWVpZdsWKF14+7ePFiLVy4MKAW1/rHP/6hN998U3/5y180cuRIn59v5cqVeuCBBzR//nwlJSUpPT1daWlpqly58ll/9vDhw1q7dq0+++wzValSRY8//rgKCgp8nhkAENyMMZ9ba7NOtY0/W/2gUaNGGjhwoG688UZNnDjRJ0WorKxML774olq3bq1OnTrpxx9/1MiRIzV06FC1adPGo6IhHVsvpFWrVrrmmmvUokUL9e/fX1deeaUOHTrk9cwAgPDAlQ0/2rFjh15//XU1bNhQb775pmrXrn3ex/zll1/08MMP65lnntHRo0fVrl07ZWZmqkKFCl5IfOxWzJw5cxQdHa0FCxYoOTnZK8cFAIQWrmwEiBo1amjkyJE6cuSImjRpogcffFBHjx4t17G++uorDR06VHXr1tWrr76qbt266dprr1VWVpbXioYkxcfHa8iQIapevbratGmjr776ymvHBgCEB65sOLJlyxa99957OnTokEaPHq1x48YpKSnpjD+zdu1avfzyy5o9e7a2bt2qjIwMtW3bVomJiX7JvGTJEi1dulTvv/++2rZt65dzAgCCw5mubFA2HLLWauPGjfr8889VXFys5ORkNW/eXA0aNFB8fLwOHjyobdu2acOGDSouLtb+/fvVuHFjNW/eXI0bN1ZUlP/fTLRs2TJ9+umnWrRokTIzM/1+fgBAYDpT2eCtrw4ZY9SwYUM1bNhQhw4dUklJiX766SetXLlShw4dUlRUlGJiYpSUlKTWrVsrMTHR+VtR27Vrp6NHj6pHjx5auXLlWa/GAABA2QgQFSpU+K14BLr27dvr559/VteuXbVy5UpVqlTJdSQAQABjgijOmTFGeXl5ioiIUO/evYPilhUAwB3KBsolIiJCffv21VdffaV7773XdRwAQACjbKDcYmJiNGDAAD388MNasGCB6zgAgABF2cB5qVWrlvLz8zV48GD9/PPPruMAAAIQZQPnLTMzU0lJSRo0aJDrKACAAETZgFfk5+frs88+07Rp01xHAQAEGMoGvKJixYrq06ePbrnlFm3ZssV1HABAAKFswGsaNmyoFi1acDsFAPAvKBvwqtzcXH399dfcTgEA/MajsmGMyTfGfGOMKTbG3HGK7VWNMe8YY74yxqw1xozwflQEgwoVKqhXr16aOHGidu7c6ToOACAAnLVsGGMiJU2RVCApXdIgY0z6SbtdJ+lra21LSV0k/bcxxnufc46g0rhxY6WkpGjkyJGuowAAAoAnVzbaSSq21n5rrT0k6VVJfU7ax0qKM8YYSVUk7ZR0xKtJEVR69Oih999/X/Pnz3cdBQDgmCdlI1nSiW8vKDn+3ImekNRc0lZJqyVNsNaWnXwgY8wYY8wKY8yK0tLSckZGMKhcubK6d++uMWPG6PDhw67jAAAc8qRsmFM8d/Inb+VJ+lJSkqRWkp4wxsT/2w9Z+4y1Nstam5WYmHiOURFsWrZsqYiICN19992uowAAHPKkbJRISjnhcT0du4JxohGS3rTHFEvaKKmZdyIiWBljVFBQoCeeeEKbNm1yHQcA4IgnZWO5pCbGmAbHJ30OlPT2Sft8J+lSSTLG1JaUJulbbwZFcKpTp44yMjI0atQo11EAAI6ctWxYa49IGi/pPUnrJM2y1q41xow1xow9vtt9kjoYY1ZL+j9JE621230VGsGlc+fOWrp0qd555x3XUQAADkR5spO1dp6keSc99/QJ32+V1MO70RAqKlasqG7dumn8+PEqKChQVJRHv3YAgBDBCqLwi8zMTBljdM8997iOAgDwM8oG/CIiIkIFBQV67LHH+KA2AAgzlA34TZ06dXThhRcyWRQAwgxlA37VpUsXffrpp5o3b97ZdwYAhATKBvzq18mi1157rY4cYUV7AAgHlA343a8rizJZFADCA2UDfvfryqJMFgWA8EDZgBO1a9dWq1atNGLECNdRAAA+RtmAM506ddLy5cs1Z84c11EAAD5E2YAzsbGx6t69u6677jodOnTIdRwAgI9QNuDUhRdeqJiYGN1xxx2uowAAfISyAad+nSw6depUrV+/3nUcAIAPUDbgXGJiorKzszVo0CBZa13HAQB4GWUDAaFTp04qKSnRI4884joKAMDLKBsICFFRUerTp4/uvfdebd682XUcAIAXUTYQMJKTk9WqVSsNHDjQdRQAgBdRNhBQOnfurOLiYj366KOuowAAvISygYASHR2toqIi3X333Vq7dq3rOAAAL6BsIOAkJycrJydHffv2ZbEvAAgBlA0EpPbt28sYo5EjR7qOAgA4T5QNBKSIiAj16dNHc+bM0XPPPec6DgDgPFA2ELCqVKmi/v3764YbbtCSJUtcxwEAlBNlAwGtfv366tatm3r37q3vv//edRwAQDlEuQ4AnE3r1q21fft2de7cWStWrFC1atVcR5Ik7dy5U8uXL9c333yj3bt3yxij2rVrKzMzU23btlVkZKTriAAQECgbCAqXXnqp5s6dq5ycHC1dulRVqlTxe4aysjK9++67mj59upYtW6Yff/xRCQkJql69umJiYmSt1b59+1RaWqrDhw+rdevWGjFihK688kpFRfG/GoDwZVx98FVWVpZdsWKF14+7ePFiLVy4kA/0CkFlZWV66623ZIzRokWL/HaF48cff9Sjjz6qV155RXv27FHLli3VsGFDJScnn7ZE7N69W+vWrdOXX34pa61uu+023XjjjYqI4M4lgNBkjPncWpt1ym2UDQSTo0ePau7cudq1a5c+/PBDXXDBBT45j7VWH3/8sSZPnqxFixapbt26ysrKUtOmTc+pMFhrtWHDBn3wwQeqUaOGZs6cqczMTJ9kBgCXzlQ2+DMLQSUyMlK9e/dW/fr11apVK7399ttePf7hw4f15JNPKiMjQ/n5+dq1a5euvvpqDR06VM2aNTvnKxPGGDVt2lSjR49WnTp11KFDBz311FNezQwAgY4byQg6xhh17dpVSUlJGjx4sIqKijRlyhTFx8eX+5ibN2/WQw89pNdee00VK1ZUVlaW+vbtq+joaK9kjoyMVKdOnZSamqo77rhDK1eu1NSpU7mtAiAsUDYQtJo1a6akpCS9//77Sk1N1fjx4zVx4kRVrlzZo58/cOCAZs6cqalTp2rVqlVq3ry5+vXrp+TkZJ9lTklJ0ejRo/XSSy9px44dmj17NoUDQMhjzgZCwqZNm/TZZ5/p+++/V4cOHVRYWKhLLrlELVq0UGRkpKy12rlzp5YuXaqPP/5YH374ob744gslJiYqMzNTGRkZio2N9Vveffv26ZVXXlHLli01Z84cCgeAoMcEUYSNnTt36ptvvtHGjRv1448/as+ePYqOjlZZWZmMMUpMTFTdunWVlJSkpk2bKi4uzlnWAwcOaMaMGcrLy9OMGTOc5QAAbzhT2eA2CkJKQkKCLr74Yl188cWSjr1d9vDhw4qIiPDa/AtviY2N1ZAhQzR9+nTdeeedevDBB11HAgCf4NotQlpERIRiYmICrmj8Ki4uTkOGDNETTzyhmTNnuo4DAD5B2QAcq1GjhoqKijR27FitXr3adRwA8DrKBhAAGjVqpPbt2+uyyy7T3r17XccBAK+ibAABokOHDqpcubKGDh3qOgoAeBVlAwgQxhgVFhbqgw8+0HPPPec6DgB4DWUDCCCVKlVSnz59dNNNN2nTpk2u4wCAV1A2gADTqFEjZWRk6PLLL1dZWZnrOABw3igbQADKzc3Vd999p//+7/92HQUAzhtlAwhAUVFR6t27t+677z5t3rzZdRwAOC+UDSBApaSkqEWLFho2bJjrKABwXigbQADLzc3VmjVrNG3aNNdRAKDcKBtAAIuJiVGvXr10++23a9euXa7jAEC5UDaAANe4cWOlpKRo7NixrqMAQLlQNoAg0L17d82dO1effPKJ6ygAcM4oG0AQiIuLU6dOnTRq1CjW3gAQdDwqG8aYfGPMN8aYYmPMHafZp4sx5ktjzFpjzCLvxgTQrl077dq1Sw8//LDrKABwTs5aNowxkZKmSCqQlC5pkDEm/aR9qkl6UlJva20LSf29HxUIb5GRkerZs6cmTZqk0tJS13EAwGOeXNloJ6nYWvuttfaQpFcl9Tlpn8GS3rTWfidJ1tpt3o0JQJIuuOACNWjQgMmiAIKKJ2UjWdKWEx6XHH/uRE0lVTfGfGiM+dwYM/xUBzLGjDHGrDDGrOAvM6B8unXrpvnz5+ujjz5yHQUAPOJJ2TCneM6e9DhKUltJvSTlSbrbGNP0337I2mestVnW2qzExMRzDgtAqlKlijp37qzRo0czWRRAUPCkbJRISjnhcT1JW0+xz3xr7V5r7XZJH0lq6Z2IAE6WlZWlPXv28EFtAIKCJ2VjuaQmxpgGxpgKkgZKevukfeZI6mSMiTLGVJJ0kaR13o0K4FeRkZEqKCjQ/fffr+3bt7uOAwBndNayYa09Imm8pPd0rEDMstauNcaMNcaMPb7POknzJa2StEzSNGvtGt/FBlC/fn01btyYyaIAAl6UJztZa+dJmnfSc0+f9PhPkv7kvWgAzqZr16566qmn9PHHH6tjx46u4wDAKbGCKBDEqlSpoi5durCyKICARtkAglzbtm21e/duJosCCFiUDSDIRUREqGfPnrr//vu1Y8cO13EA4N9QNoAQkJKSosaNG2vMmDGuowDAv6FsACHi0ksv1XvvvadFi/gcRACBhbIBhIjKlSv/Nln06NGjruMAwG8oG0AIadu2rQ4ePKi77rrLdRQA+A1lAwghERERKiws1GOPPaZ161jEF0BgoGwAIaZWrVpq3769Bg4cyNobAAICZQMIQTk5Odq2bZseeugh11EAgLIBhKLIyEj17t1bDzzwgNavX+86DoAwR9kAQlRSUpKys7NVVFTEu1MAOEXZAEJYp06dtHfvXk2YMMF1FABhjLIBhLDIyEgVFRVpxowZmj9/vus4AMIUZQMIcdWqVVPPnj01dOhQ/fDDD67jAAhDlA0gDLRo0UJNmzZVXl6ejhw54joOgDBD2QDCRLdu3bRv3z4NHjzYdRQAYYayAYSJyMhI9evXT//3f/+nO++803UcAGEkynUAAP5TqVIlDR06VE8++aRq1KihW2+91Wme4uJivfnmm1q8eLE2bNig7du36+DBg5KkuLg41a5dW9nZ2erXr5+6d++uiAj+PgKCEWUDCDMJCQkaMmSI7r33XlWsWFHXXXedX8+/efNmPf7445ozZ45KSkrUoEEDJScnKzs7W1WrVlVsbKwk6ZdfftHOnTu1atUqvfHGG4qJidGIESN01113qWLFin7NDOD8GGutkxNnZWXZFStWeP24ixcv1sKFC+XqvwsIFlu3btXMmTN122236Z577vHpuay1evfdd/XHP/5RS5cuVdOmTdWiRQs1atRIUVFn/5unrKxMGzdu1CeffKI9e/Zo0qRJGjVqlE8zAzg3xpjPrbVZp9xG2QDC1/bt2/XSSy+pZ8+eev755z164T8XBw8e1OOPP66nnnpKO3bsUFZWltq0aaPKlSuX63jWWm3YsEHz5s1Tdna2XnvtNVWvXt2rmQGUz5nKBjdAgTBWs2ZNjRo1SkuXLlVGRobWrl3rleNu3bpV48aNU+3atTVlyhRlZWXphhtuUKdOncpdNCTJGKOmTZtq7Nix2rZtmzIyMrRmzRqvZAbgO5QNIMxVqVJFQ4YMUXJystq1a6fx48drz54953wca63mz5+vvLw8NWrUSEuWLNHAgQN11VVXKT09XZGRkV7LHBsbq759+yo9PV05OTlasGCB144NwPuYIApAERER6tSpk9LT07VgwQIlJSWpf//+GjdunLKyTnlVVJJ05MgRLViwQK+//rrmzp2rI0eOqGXLlrruuusUFxfn08zGGOXk5KhatWoqKirS7NmzlZeX59NzAigfygaA39SoUUNXXHGFfvrpJ33xxRfq2rWrKlSooEaNGiklJUXx8fE6ePCgduzYoY0bN2rLli2Ki4tTo0aNVFhYqHr16skY49fMLVq0UEREhPr166e5c+eqS5cufj0/gLNjgiiA0yorK9OOHTv0ww8/aM+ePTp48KAiIyNVsWJF1ahRQ7Vr11aVKlVcx5QkrVq1Sh988IE+/fRTNW/e3HUcIOycaYIoVzYAnFZERIQSExOVmJjoOspZZWZm6p///Ke6deumL774QrVq1XIdCcBxTBAFEDI6duyopKQkFRQUqKyszHUcAMdRNgCEDGOM8vPztW3bNl177bWu4wA4jrIBIKRERUWpf//+mjlzpl555RXXcQCIsgEgBFWtWlV9+/bVtddeqy1btriOA4Q9ygaAkNS4cWOlp6fr8ssvZ/4G4BhlA0DI6tq1qzZt2qRJkya5jgKENcoGgJAVHR2toqIiTZ48WatXr3YdBwhblA0AIa1u3bq66KKLNGTIEG6nAI5QNgCEvJycHG3btk0PPfSQ6yhAWKJsAAh5kZGRKiws1IMPPqjNmze7jgOEHcoGgLBQr149tWjRQsOHD3cdBQg7lA0AYSM3N1erVq3SCy+84DoKEFYoGwDCRkxMjAoKCnTLLbdo7969ruMAYYOyASCspKWlqWbNmpowYYLrKEDYoGwACDt5eXl65ZVXWHsD8BPKBoCwU716dbVr105XX3216yhAWKBsAAhLOTk5+vbbbzV9+nTXUYCQR9kAEJaio6OVl5eniRMnMlkU8DHKBoCwlZaWpoSEBN10002uowAhzaOyYYzJN8Z8Y4wpNsbccYb9so0xR40xl3svIgD4Tl5enl5++WWtXbvWdRQgZJ21bBhjIiVNkVQgKV3SIGNM+mn2e0jSe94OCQC+kpCQoOzsbCaLAj7kyZWNdpKKrbXfWmsPSXpVUp9T7He9pDckbfNiPgDwuY4dO2rDhg16/vnnXUcBQpInZSNZ0pYTHpccf+43xphkSUWSnvZeNADwj+joaOXn5+u2227Tvn37XMcBQo4nZcOc4jl70uM/S5porT16xgMZM8YYs8IYs6K0tNTDiADge2lpaapRo4ZuvPFG11GAkONJ2SiRlHLC43qStp60T5akV40xmyRdLulJY0zfkw9krX3GWptlrc1KTEwsX2IA8JFfJ4uuWbPGdRQgpHhSNpZLamKMaWCMqSBpoKS3T9zBWtvAWptqrU2VNFvSOGvtX70dFgB8qXr16mrfvr1GjBjhOgoQUs5aNqy1RySN17F3mayTNMtau9YYM9YYM9bXAQHAnzp06KCNGzfq2WefdR0FCBlRnuxkrZ0nad5Jz51yMqi19qrzjwUAbkRFRamgoEC33367BgwYoLi4ONeRgKDHCqIAcJLGjRurbt26uv76611HAUICZQMATqFHjx56/fXXtXTpUtdRgKBH2QCAU6hatao6d+6soUOH6ujRM76rH8BZUDYA4DSys7N1+PBh3Xnnna6jAEGNsgEApxEREaHCwkJNmTJFq1atch0HCFqUDQA4g5o1ayonJ0cDBw5UWVmZ6zhAUKJsAMBZtG/fXrt37+Z2ClBOlA0AOIvIyEgVFRXpiSee0EcffeQ6DhB0KBsA4IGaNWuqR48eGjhwoHbv3u06DhBUKBsA4KGWLVsqMTFRAwYMcB0FCCqUDQDwkDFGPXv21PLly/XAAw+4jgMEDcoGAJyD2NhYDRw4UJMmTdLcuXNdxwGCAmUDAM5RrVq11KdPHw0ZMkTr1q1zHQcIeB596isA4F+lpaVp586dys3N1bJly1S/fn3XkWSt1bp16/TJJ59o48aN+uc//yljjBITE5Wenq5LL71UCQkJrmMiDFE2AKCcLr74Yh04cEAdO3bU8uXLVbt2bb9n+OWXX/Tiiy/qtdde05dffqnDhw8rKSlJ8fHxio2NlbVW+/fv17PPPquffvpJDRo00ODBg3XTTTcpLi7O73kRnigbAHAeunTpooMHD6pdu3ZavHix365wLFq0SJMnT9aHH36oOnXqqFmzZho+fLgSEhJkjDnlzxw8eFDffvutXnrpJf3pT3/S2LFjdf/99ysmJsYvmRG+jLXWyYmzsrLsihUrvH7cxYsXa+HChXL13wUg/Fhr9eGHH+rrr7/W/Pnz1aZNG5+c59ChQ3ryySf11FNPaevWrWrbtq3atGmjqlWrnvOxfvjhB73//vs6cOCAXnzxRXXt2tUHiRFOjDGfW2uzTrWNCaIAcJ6MMcrNzVW7du3UuXNnPf300149/k8//aTrr79etWvX1qOPPqqMjAzdeOONys3NLVfRkKS6detq2LBhys7OVmFhoSZOnMgfafAZbqMAgJdkZ2erTp06uvPOOzV79mxNnz5dKSkp5T7eihUr9F//9V9asGCBmjZtqv79+ys5OdlreY0xat26terXr6/p06dr3bp1euONNxQdHe21cwASVzYAwKtSUlI0duxY7d+/X82bN9fVV1+tkpISj3++tLRUf/jDH5SWlqbOnTvr559/1tixY1VUVOTVonGiGjVq6Oqrr9batWt1ySWX6ODBgz45D8IXczYAwEe2b9+uJUuWaPXq1WrevLny8/PVoUMHtWjRQjVr1tSBAwdUUlKilStXasmSJVq8eLE2bdqkRo0aKTMzU02aNFFUlP8uQB8+fFizZs1SjRo1tGjRIlWoUMFv50bwO9OcDcoGAPjYgQMHtH79em3atEk//vijfv75Zx06dEhRUVGqWLGiatSooVq1aik1NVWpqalOX+QPHz6s1157TRdccIH+9re/KSKCC+DwzJnKBnM2AMDHYmNjlZmZqczMTNdRzio6Olr9+/fXjBkzNGbMGE2bNs11JIQAKisA4F/ExMRo0KBBmj17tv70pz+5joMQQNkAAPyb+Ph4DRw4UPfee68WLVrkOg6CHGUDAHBKdevWVffu3TVgwACVlpa6joMgRtkAAJxW69atlZKSor59+6qsrMx1HAQpygYA4Izy8/NVXFysSZMmuY6CIEXZAACcUXR0tPr27avJkydr9erVruMgCFE2AABnlZSUpIsuukgDBgzQ0aNHXcdBkKFsAAA8kpOTo7179+q2225zHQVBhrIBAPBIZGSk+vTpo6lTp2rVqlWu4yCIUDYAAB6rWbOm2rdvr2HDhvHuFHiMsgEAOCcdOnTQTz/9xOqi8BhlAwBwTqKionTZZZfpgQceUElJies4CAKUDQDAOUtJSVF6erquuuoq11EQBCgbAIByyc3N1bJly/TWW2+5joIAR9kAAJRLbGysunXrphtuuEGHDx92HQcBjLIBACi3zMxMRUVF6T/+4z9cR0EAo2wAAMrNGKP8/Hw99dRT2rRpk+s4CFCUDQDAealTp44yMjI0cuRI11EQoCgbAIDz1qVLFy1fvlzvvPOO6ygIQJQNAMB5i42N1aWXXqrx48czWRT/hrIBAPCKli1byhij//zP/3QdBQGGsgEA8ApjjAoKCvTYY4/pu+++cx0HAYSyAQDwmjp16igzM5PJovgXlA0AgFd17txZS5cu1dy5c11HQYCgbAAAvCo2Nlbdu3fXuHHjdOTIEddxEAA8KhvGmHxjzDfGmGJjzB2n2D7EGLPq+NenxpiW3o8KAAgWGRkZio6O1u9//3vXURAAzlo2jDGRkqZIKpCULmmQMSb9pN02Supsrc2UdJ+kZ7wdFAAQPH5dWXTKlCmsLAqPrmy0k1Rsrf3WWntI0quS+py4g7X2U2vtruMPl0iq592YAIBgU6tWLbVu3VojRoxwHQWOeVI2kiVtOeFxyfHnTmekpHdPtcEYM8YYs8IYs6K0tNTzlACAoNSpUyd98cUXmj17tusocMiTsmFO8Zw95Y7G5OpY2Zh4qu3W2mestVnW2qzExETPUwIAglJMTIx69Oih66+/XgcOHHAdB454UjZKJKWc8LiepK0n72SMyZQ0TVIfa+0O78QDAAS79PR0xcXFacKECa6jwBFPysZySU2MMQ2MMRUkDZT09ok7GGPqS3pT0jBr7XrvxwQABCtjjHr16qWXXnpJn332mes4cOCsZcNae0TSeEnvSVonaZa1dq0xZqwxZuzx3e6RVEPSk8aYL40xK3yWGAAQdKpVq6bc3FwNGTKEtTfCkEfrbFhr51lrm1prG1lrHzj+3NPW2qePfz/KWlvdWtvq+FeWL0MDAIJP27ZtJUk333yz4yTwN1YQBQD4RUREhAoLC/Xss89q2bJlruPAjygbAAC/SUhIUJcuXTRgwADenRJGKBsAAL/Kzs5WbGysrrzyStdR4CeUDQCAXxlj1Lt3b82fP1/PP/+86zjwA8oGAMDvKleurKKiIl1//fVav54VE0IdZQMA4ESDBg3Url079ejRQ7/88ovrOPAhygYAwJmOHTsqPj5eBQUFKisrcx0HPkLZAAA48+v8jY0bN+rqq692Hee0jh49ymJk54GyAQBwKjo6WgMHDtTcuXN10003Oc1SVlamjz76SDfffLPat2+vpKQkxcbGKioqShUqVFCVKlWUmpqqwsJCTZs2TQcPHnSaN1hEuQ4AAEBcXJyGDRumGTNmKDY2Vg8++KBfz79582Y9+uijmjVrlvbv369GjRqpfv36atWqlapXr64KFSrIWqv9+/dr+/bt2rJli/7whz/opptuUv/+/TV58mTVqlXLr5mDCWUDABAQqlWrpuHDh+vpp59WaWmpnnnmGUVE+O4CfFlZmd5991098sgj+uyzz5Samqpu3bqpQYMGpz1vlSpVfru60alTJ5WWlmrx4sVq2LChbrvtNt19990+zRysjLXWyYmzsrLsihXe/7y2xYsXa+HChXL13wUAOD979uzRzJkz1bx5c7311luKj4/36vH37t2rxx57TM8++6y2bdum1q1bq23btqpatWq5j7l161a9/fbbqlWrlt555x3Vr1/fi4mDgzHm89N9Nhr1CwAQUOLi4nTVVVdp+/btatasmRYuXOiV4/7973/XsGHDVLduXT3zzDNq3bq1JkyYoK5du55X0ZCkpKQkjR49WlWqVFHr1q21YMECr2QOFZQNAEDAiYmJUVFRkbKzs1VYWKi+fftqy5Yt53ycQ4cO6X/+53+UlZWlNm3aaP369Ro6dKiuuuoqXXjhhYqK8t5sgsjISPXo0UO5ubnq27evXn31Va8dO9gxZwMAEJCMMWrTpo3S0tL04YcfqmnTpurcubNGjhypvn37Kjo6+pQ/t3PnTv31r3/V7Nmz9fHHH6tatWpq1aqVevTooZiYGJ/nzsjIUOXKlTVq1Cjt379fI0aM8Pk5Ax1zNgAAQWHPnj364osv9Pe//107d+5UUlKS6tWrp0qVKunIkSPauXOnSkpKtGvXLiUnJ6tJkyZKS0tTzZo1neT97rvv9Nprr+nFF19UUVGRkwz+dKY5G5QNAEDQ2bt3r0pLS7Vr1y4dPnxYERERqlSpkhISEpSYmKjIyEjXESVJGzZs0Jw5c/S///u/uuSSS1zH8akzlQ1uowAAgk7lypVVuXJlpaamuo5yRk2aNFH37t3Vp08frVy5Ug0aNHAdyQkmiAIA4EMtW7ZUZmam8vPzw3bFUcoGAAA+1qVLF1lrdcUVV7iO4gRlAwAAH4uIiNDvfvc7LV68WH/5y19cx/E7ygYAAH5QsWJFFRUV6fe//72+/vpr13H8irIBAICfpKSkqF27durXr19YfWQ9ZQMAAD/q2LGj9u/frxtuuMF1FL+hbAAA4EeRkZHq27evnn/+eX3yySeu4/gFZQMAAD+rXr26OnfurOHDh4fF7RTKBgAADmRnZ+vIkSO6/fbbXUfxOcoGAAAOREREqLCwUFOnTtXq1atdx/EpygYAAI4kJiaqffv2Gjx4sMrKylzH8RnKBgAADuXk5Ki0tFSPPvqo6yg+Q9kAAMChyMhIFRQU6L777tP27dtdx/EJygYAAI6lpqYqNTVV11xzjesoPkHZAAAgAHTr1k3z58/XokWLXEfxOsoGAAABIC4uTpdcconGjBkTcpNFKRsAAASI7Oxs7d69Ww899JDrKF5F2QAAIEBERkaqZ8+emjx5srZt2+Y6jtdQNgAACCD169dX48aNQ2qyKGUDAIAA07VrV73//vv66KOPXEfxCsoGAAABpkqVKurSpYtGjx4dEpNFKRsAAASgtm3bas+ePXr44YddRzlvlA0AAAJQRESEevbsqfvvv1+lpaWu45wXygYAAAEqJSVFTZs2DfrJopQNAAACWNeuXfW3v/0tqCeLUjYAAAhglStXVm5urkaNGhW0k0UpGwAABLg2bdpo7969+uMf/+g6SrlQNgAACHC/ThadNGmSfvrpJ9dxzhllAwCAIFCvXj2lpaVp9OjRrqOcM4/KhjEm3xjzjTGm2Bhzxym2G2PMY8e3rzLGtPF+VAAAwlvXrl21cOFCzZs3z3WUc3LWsmGMiZQ0RVKBpHRJg4wx6SftViCpyfGvMZKe8nJOAADCXqVKldSjRw+NGjVK+/fvdx3HY55c2Wgnqdha+6219pCkVyX1OWmfPpJesMcskVTNGFPXy1kBAAh7GRkZio+P17hx41xH8ZgnZSNZ0pYTHpccf+5c9wEAAOfJGKNevXpp1qxZWrRokes4HonyYB9ziudsOfaRMWaMjt1mkaRfjDHfeHD+8qgpabuPjo1/xVj7F+PtP4y1/zDW5dSlS5fy/JivxvuC023wpGyUSEo54XE9SVvLsY+stc9IesaDc54XY8wKa22Wr88DxtrfGG//Yaz9h7H2Lxfj7cltlOWSmhhjGhhjKkgaKOntk/Z5W9Lw4+9KaS/pn9baH7ycFQAABKGzXtmw1h4xxoyX9J6kSEnPWWvXGmPGHt/+tKR5knpKKpa0T9II30UGAADBxJPbKLLWztOxQnHic0+f8L2VdJ13o50Xn9+qwW8Ya/9ivP2HsfYfxtq//D7e5lhPAAAA8A2WKwcAAD4V1GWDZdT9x4OxHnJ8jFcZYz41xrR0kTMUnG2sT9gv2xhz1BhzuT/zhRpPxtsY08UY86UxZq0xJjgWNghAHvw7UtUY844x5qvjY838v3IyxjxnjNlmjFlzmu3+fX201gbll45NVv2HpIaSKkj6SlL6Sfv0lPSujq0D0l7SUte5g/HLw7HuIKn68e8LGGvfjfUJ+32gY3OpLnedO1i/PPzdribpa0n1jz+u5Tp3MH55ONb/Iemh498nStopqYLr7MH4JekSSW0krTnNdr++PgbzlQ2WUfefs461tfZTa+2u4w+X6NhaKzh3nvxeS9L1kt6QtM2f4UKQJ+M9WNKb1trvJMlay5iXjydjbSXFGWOMpCo6VjaO+DdmaLDWfqRj43c6fn19DOaywTLq/nOu4zhSxxozzt1Zx9oYkyypSNLTwvny5He7qaTqxpgPjTGfG2OG+y1daPFkrJ+Q1FzHFoVcLWmCtbbMP/HCjl9fHz1662uA8toy6jgrj8fRGJOrY2Wjo08ThS5PxvrPkiZaa48e+wMQ58GT8Y6S1FbSpZIqSvrMGLPEWrve1+FCjCdjnSfpS0ldJTWS9L4xZrG1drePs4Ujv74+BnPZ8Noy6jgrj8bRGJMpaZqkAmvtDj9lCzWejHWWpFePF42aknoaY45Ya//ql4ShxdN/R7Zba/dK2muM+UhSS0mUjXPjyViPkDTZHptUUGyM2SipmaRl/okYVvz6+hjMt1FYRt1/zjrWxpj6kt6UNIy/+M7LWcfaWtvAWptqrU2VNFvSOIpGuXny78gcSZ2MMVHGmEqSLpK0zs85Q4EnY/2djl1BkjGmtqQ0Sd/6NWX48OvrY9Be2bAso+43Ho71PZJqSHry+F/cRywfrHTOPBxreIkn422tXWeMmS9plaQySdOstad8OyFOz8Pf7fskzTDGrNaxy/wTrbV8Gmw5GGNekdRFUk1jTImk/5QULbl5fWQFUQAA4FPBfBsFAAAEAcoGAADwKcoGAADwKcoGAADwKcoGAADwKcoGAADwKcoGAADwKcoGAADwqf8HVxyw40F5AnAAAAAASUVORK5CYII=\n", 58 "text/plain": [ 59 "<Figure size 648x432 with 1 Axes>" 60 ] 61 }, 62 "metadata": { 63 "needs_background": "light" 64 }, 65 "output_type": "display_data" 66 } 67 ], 68 "source": [ 69 "from matplotlib.pylab import plt\n", 70 "import numpy as np\n", 71 "plt.rcParams[\"figure.figsize\"] = (9,6)\n", 72 "\n", 73 "fig = plt.figure()\n", 74 "x_grid = np.linspace(0, 1., 1000)\n", 75 "plt.plot(x_grid, (1 - x_grid + 0.05 * np.cos(11 * np.pi * x_grid)), 'k-', linewidth=1)\n", 76 "plt.fill_between(x_grid, -1., (1 - x_grid + 0.05 * np.cos(11 * np.pi * x_grid)), color='gray')\n", 77 "plt.ylim(0, None);" 78 ] 79 }, 80 { 81 "cell_type": "markdown", 82 "id": "a02d74cf", 83 "metadata": {}, 84 "source": [ 85 "We want the particle to bounce off both the wavy ramp and the ground (i.e., $x$ axis). Thus, we will have to define two different event equations and two different callback functions to implement the bouncing behaviour. Let us begin with the event equations:" 86 ] 87 }, 88 { 89 "cell_type": "code", 90 "execution_count": 3, 91 "id": "04c9f0e8", 92 "metadata": {}, 93 "outputs": [], 94 "source": [ 95 "# Event equation for the ramp.\n", 96 "eq_w_curve = y - (1. - x + 0.05 * hy.cos(11 * np.pi * x))\n", 97 "\n", 98 "# Event equation for the ground.\n", 99 "eq_bottom = y" 100 ] 101 }, 102 { 103 "cell_type": "markdown", 104 "id": "4479fba2", 105 "metadata": {}, 106 "source": [ 107 "Next, we need to implement the callbacks. Because we are dealing with inelastic collisions, when the particle hits a surface it will lose part of its kinetic energy. The [coefficient of restitution](https://en.wikipedia.org/wiki/Coefficient_of_restitution) $C_R \\in \\left[0, 1\\right]$ regulates how much energy is lost at every bounce.\n", 108 "\n", 109 "Additionally, because of the loss of energy, the particle will experience increasingly frequent collisions if it gets stuck in one of the local minima of the curve or when it starts bouncing off the ground. In practice, we want to stop the simulation once the collisions become too frequent.\n", 110 "\n", 111 "Let us see the code:" 112 ] 113 }, 114 { 115 "cell_type": "code", 116 "execution_count": 4, 117 "id": "c6a7ab81", 118 "metadata": {}, 119 "outputs": [], 120 "source": [ 121 "# The coefficient of restitution.\n", 122 "CR = .8\n", 123 "\n", 124 "# Global variables to track\n", 125 "# the time of the last collision\n", 126 "# and the collision points.\n", 127 "last_t = 0.\n", 128 "bounce_points = []\n", 129 "\n", 130 "# Callback for bouncing against the curve.\n", 131 "def cb_w_curve(ta, mr, d_sgn):\n", 132 " global last_t\n", 133 "\n", 134 " # If the last collision happened\n", 135 " # too recently, return False to\n", 136 " # stop the simulation.\n", 137 " if ta.time - last_t < 1e-6:\n", 138 " return False\n", 139 "\n", 140 " # Update last_t.\n", 141 " last_t = ta.time\n", 142 "\n", 143 " # Update bounce_points.\n", 144 " x, y = ta.state[0:2]\n", 145 " bounce_points.append((x, y))\n", 146 " \n", 147 " # Compute the normal unit vector\n", 148 " # using the gradient of the event\n", 149 " # equation.\n", 150 " grad = np.array([1+0.05*11*np.pi*np.sin(11*np.pi*x), 1])\n", 151 " grad_uvec = grad / np.linalg.norm(grad)\n", 152 " \n", 153 " # Compute the component of the velocity\n", 154 " # across the normal vector.\n", 155 " xy_vel = ta.state[2:4]\n", 156 " vproj = np.dot(xy_vel, grad_uvec)\n", 157 " \n", 158 " # Flip it and rescale it according\n", 159 " # to the coefficient of restitution.\n", 160 " Dv = -vproj*grad_uvec\n", 161 " xy_vel += (1. + CR) * Dv\n", 162 "\n", 163 " return True\n", 164 "\n", 165 "# Callback for bouncing off the ground.\n", 166 "def cb_bottom(ta, mr, d_sgn):\n", 167 " global last_t\n", 168 " \n", 169 " # If the last collision happened\n", 170 " # too recently, return False to\n", 171 " # stop the simulation.\n", 172 " if ta.time - last_t < 1e-6:\n", 173 " return False\n", 174 "\n", 175 " # Update last_t.\n", 176 " last_t = ta.time\n", 177 "\n", 178 " # Update bounce_points.\n", 179 " x, y = ta.state[0:2]\n", 180 " bounce_points.append((x, y))\n", 181 " \n", 182 " # Flip the y component of the velocity,\n", 183 " # and rescale it according\n", 184 " # to the coefficient of restitution.\n", 185 " ta.state[3] = -CR*ta.state[3]\n", 186 " \n", 187 " return True" 188 ] 189 }, 190 { 191 "cell_type": "markdown", 192 "id": "ba77fb48", 193 "metadata": {}, 194 "source": [ 195 "We are now ready to create the integrator object. As initial conditions, we will put the particle at rest above the ramp on the $y$ axis:" 196 ] 197 }, 198 { 199 "cell_type": "code", 200 "execution_count": 5, 201 "id": "285b4776", 202 "metadata": {}, 203 "outputs": [], 204 "source": [ 205 "ta = hy.taylor_adaptive(eqns, [0, 1.2, 0, 0], t_events = [hy.t_event(eq_w_curve, callback = cb_w_curve, direction = hy.event_direction.negative),\n", 206 " hy.t_event(eq_bottom, callback = cb_bottom, direction = hy.event_direction.negative)])" 207 ] 208 }, 209 { 210 "cell_type": "markdown", 211 "id": "e7b95f4e", 212 "metadata": {}, 213 "source": [ 214 "Note how, as explained in the [Keplerian billiard](<./The Keplerian billiard.ipynb>) example, we assigned specific directions to the collision events. This is done in order to avoid triggering spurious double-bounce events that would lead to the particle penetrating through the ramp or the ground.\n", 215 "\n", 216 "We can now integrate the system for a few time units:" 217 ] 218 }, 219 { 220 "cell_type": "code", 221 "execution_count": 6, 222 "id": "a1fab993", 223 "metadata": {}, 224 "outputs": [], 225 "source": [ 226 "t_grid = np.linspace(0, 10, 10000)\n", 227 "oc, _, _, _, res = ta.propagate_grid(t_grid)\n", 228 "\n", 229 "# Transform bounce_points in a NumPy array\n", 230 "# for ease of use.\n", 231 "b_points = np.array(bounce_points)" 232 ] 233 }, 234 { 235 "cell_type": "markdown", 236 "id": "1c503814", 237 "metadata": {}, 238 "source": [ 239 "Let us take a look at the trajectory:" 240 ] 241 }, 242 { 243 "cell_type": "code", 244 "execution_count": 7, 245 "id": "cff0eff9", 246 "metadata": {}, 247 "outputs": [ 248 { 249 "data": { 250 "image/png": "iVBORw0KGgoAAAANSUhEUgAAAsMAAAIICAYAAACLo+M1AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAB41UlEQVR4nO3dd3gc1dk28Pts06qXberNknvvxgVsY4wBx/QYQglJ6CXkSwLhTQK8SUh9U0iAhB4goRdjwIDBGNyLZFuyLNmWbEtW771tm++PlWSVlbSSZou19++6uOLdmTlztFnM7eNnniMkSQIRERERkT9SeHsCRERERETewjBMRERERH6LYZiIiIiI/BbDMBERERH5LYZhIiIiIvJbDMNERERE5LdU3rqxXq+XkpOTvXV7IiIiIvITmZmZNZIkGZwd81oYTk5ORkZGhrduT0RERER+QghRNNgxlkkQERERkd9iGCYiIiIiv8UwTERERER+i2GYiIiIiPwWwzARERER+S2GYSIiIiLyWwzDREREROS3GIaJiIiIyG8xDBMRERGR32IYJiIiIiK/xTBMRERERH6LYZiIiIiI/BbDMBERERH5LYZhIiIiIvJbDMNERERE5LcYhomIiIjIbzEMExEREZHfYhgmIiIiIr/FMExEREREfothmIiIiIj8FsMwEREREfkthmEiIiIi8lsMw0RERETktxiGiYiIiMhvMQwTERERkd9iGCYiIiIiv8UwTERERER+i2GYiIiIiPwWwzARERER+S2GYSIiIiLyWwzDREREROS3GIaJiIiIyG8xDBMRERGR32IYJiIiIiK/xTBMRERERH5r2DAshHhJCFElhMgZ5Ph3hBDZXf/sEULMkn+aRERERETyc2Vl+N8ALh3i+BkAF0qSNBPArwE8J8O8iIiIiIjcbtgwLEnSDgB1QxzfI0lSfdfLfQDiZZqb7Kw2O8ob29Fmtnp7KkRERETkA+SuGf4+gE9lHlM2Vc2dWPK7r/BRVpm3p0JEREREPkAl10BCiJVwhOFlQ5xzB4A7ACAxMVGuWxMRERERjYosK8NCiJkAXgCwQZKk2sHOkyTpOUmS5kuSNN9gMMhxayIiIiKiURtzGBZCJAJ4H8DNkiSdHPuUiIiIiIg8Y9gyCSHEGwAuAqAXQpQAeAyAGgAkSfoXgEcB6AA8I4QAAKskSfPdNWEiIiIiIrkMG4YlSbphmOM/APAD2WZEREREROQh3IGOiIiIiPwWwzARERER+S2GYSIiIiLyWwzDREREROS3GIaJiIiIyG8xDBMRERGR32IYJiIiIiK/xTBMRERERH6LYZiIiIiI/BbDMBERERH5LYZhIiIiIvJbDMNERERE5Lf8MgxLkrdnQERERES+wK/CsBDengERERER+RK/CsNERERERL0xDBMRERGR32IYJiIiIiK/xTBMRERERH6LYZiIiIiI/BbDMBERERH5LYZhIiIiIvJbDMNERERE5LcYhomIiIjIbzEMExEREZHfYhgmIiIiIr/FMExEREREfothmIiIiIj8FsMwEREREfkthmEiIiIi8lsMw0RERETktxiGiYiIiMhv+WUYlrw9ASIiIiLyCX4VhgWEt6dARERERD7Er8IwEREREVFvDMNERERE5LcYhomIiIjIbzEMExEREZHfYhgmIiIiIr/FMExEREREfothmIiIiIj8FsMwEREREfkthmEiIiIi8lsMw0RERETktxiGiYiIiMhvMQwTERERkd9iGCYiIiIiv8UwTERERER+i2GYiIiIiPyWX4Vhbd572KV5ABu3zAT+Oh3IftvbUyIiIiIiL1J5ewIek/02wr78MSIU7Y7XjcXARw84fj3zeu/Ni4iIiIi8xn9Whrf9Cgpre9/3LO3Atl95Zz5ERERE5HX+E4YbS5y+bW8swa78Gg9PhoiIiIh8gf+E4fB4p29XQo+bXtyPm1/cj/zKZg9PioiIiIi8yX/C8OpHAXVg3/fUgdBf+Rv84vIpyC5pxOV/34WnvsqHxWb3zhyJiIiIyKP8JwzPvB5Y/3c0inDYJaDGGgis/zvUszfiB8tTse3HF2LNNBP+b+tJ3PDcPlQ2dXh7xkRERETkZv4ThgFg5vV4Luh+3FZ4Jeb/R9uni4Q+JABP3zgXT26cjdzyJlz+913Yf7rWi5MlIiIiInfzrzDcJTExEeXl5SgrKxtwbMPsOGy6dynCtCrc/OIBbDla7oUZEhEREZEn+GUYVqlUSEhIwIcffuj0+ERTKD64Zylmxofj3tcP4bV9RR6eIRERERF5gl+GYQBISEjAtm3bBj0eHqTGa99fhNWTjfjlphy8sqfQc5MjIiIiIo/w2zCclJSEzMzMIc8J1Cjxz5vmYc1UEx7bfAxvHTzrodkRERERkSf4bRiOi4tDaWkpGhoahjxPrVTgqRvn4MKJBvzs/aOsISYiIiIaR/w2DGs0GhiNxiFLJboFqJR49uZ5mJsYiR+9dQRHihvcP0EiIiIicju/DcMAEBsbi507d7p0rlatxHM3z4MxLAA/eCUDpQ3tbp4dEREREbmbX4fh6OjoYeuGe9OFBOClWxeg02LDD17JQLvZ5sbZEREREZG7+XUYjo2NRX5+/oiuSTeF4h83zsHxiiY8vvmYm2ZGRERERJ7g12HYYDCgrq4ONTU1I7ruoklG3HtRGt7KKMa7mSVumh0RERERuduwYVgI8ZIQokoIkTPIcSGE+LsQokAIkS2EmCv/NN1DpVLBZDJh+/btI772wYvTsTg1Cr/YdBQnK5vdMDsiIiIicjdXVob/DeDSIY6vA5De9c8dAP459ml5TkxMDPbs2TPi61RKBf6+cQ6CNSr86K0jMFvtbpgdEREREbnTsGFYkqQdAOqGOGUDgFclh30AIoQQMXJN0N2MRiOys7NHd22YFr+9egaOlTXhqe0FMs+MiIiIiNxNjprhOADFvV6XdL03gBDiDiFEhhAio7q6WoZbj53RaMSpU6dGff3aadG4em4cnt5egCz2HyYiIiI6r8gRhoWT9yRnJ0qS9JwkSfMlSZpvMBhkuPXYGY1GlJeXw2q1jnqMx9ZPgzE0AD9+JwudVrZbIyIiIjpfyBGGSwAk9HodD6BMhnE9IjAwEIGBgTh8+PCoxwgPVON3V89AQVULnvvmtIyzIyIiIiJ3kiMMbwZwS1dXicUAGiVJKpdhXI8xmUzYt2/fmMa4aJIRl8+MwT+2F6CwplWmmRERERGRO7nSWu0NAHsBTBJClAghvi+EuEsIcVfXKVsAnAZQAOB5APe4bbZuYjAYxrQy3O3RK6YiQKnALz/MgSQ5rRQhIiIiIh+iGu4ESZJuGOa4BOBe2WbkBUajEceOjX03OVOYFj9ZOwmPbT6Gj7LL8a1ZsTLMjoiIiIjcxa93oOtmMplQVFQky1g3LU7CjLhwPPFJLtrMo38oj4iIiIjcj2EYgF6vR21tLVpbx17rq1QIPLZ+KiqbOvHcDj5MR0REROTLGIbh2JY5MjISBw4ckGW8+clRuHxGDJ795jQqGjtkGZOIiIiI5Mcw3MVkMuHgwYOyjfezdZNhs0v4v60nZBuTiIiIiOTFMNxFr9cjKytLtvESooJw27JkvHeoBDmljbKNS0RERETyYRjuYjAYcPz4cVnHvHdlGiKDNPjDZ/KOS0RERETyYBjuYjQacfbsWVnHDNOqcc9FE7Azvwb7TtfKOjYRERERjR3DcJeoqCg0NDSgublZ1nFvWpwEY2gA/rz1BDfiICIiIvIxDMNdlEol9Hq9bB0lumnVSty/Kg0HC+uxI79G1rGJiIiIaGwYhnsxGo3IyMiQfdxvL0hEXEQgV4eJiIiIfAzDcC9yd5ToplEp8MOL05Fd0ogv86pkH5+IiIiIRodhuBd3dJTodvWcOCREBeKp7QVcHSYiIiLyEQzDvRiNRhQVFbllbJVSgbsunICs4gbsPcXOEkRERES+gGG4l8jISDQ3N8veUaLbNXPjYQwNwNNfF7hlfCIiIiIaGYbhXhQKBfR6Pfbv3++W8bVqJW5fnordBbU4fLbeLfcgIiIiItcxDPcTHR2NgwcPum38GxclIjxQjWe+PuW2exARERGRaxiG+9Hr9cjOznbb+MEBKty2NBlf5FbiRIV7yjGIiIiIyDUMw/3o9Xrk5eW59R7fvSAZgWolXtx12q33ISIiIqKhMQz3486OEt0igjS4dl48Nh0pQ01Lp1vvRURERESDYxjuJzw8HG1tbaivd+8Dbt9dmgyz1Y7/7HNv8CYiIiKiwTEM96NQKGAwGNyyLXNvEwwhWDXZiP/sK0Kn1ebWexERERGRcwzDTkRHR+PAgQNuv8/3lqagpsWMzUfK3H4vIiIiIhqIYdgJg8GArKwst99naZoOk0yheGl3IbdoJiIiIvIChmEndDodjh8/7vb7CCHwvWXJyCtvwt7T3KKZiIiIyNMYhp3wREeJbhtmxyEySM0H6YiIiIi8gGHYibCwMJjNZtTU1Lj9Xlq1EtfNT8DWY5Woaupw+/2IiIiI6ByGYSeEEDAajdi3b59H7nfDwkRY7RLeOljskfsRERERkQPD8CBMJpNHOkoAQIo+GMvT9XjjwFnY7HyQjoiIiMhTGIYHodPpkJ2d7bH7fWdRIsoaO/D1iSqP3ZOIiIjI3zEMD8JgMODEiRMeu9/qKSYYQwP4IB0RERGRBzEMD8JoNOLs2bMe6/+rViqwcWEivj5ZjeK6No/cczhFta14N7MEdpZuEBER0TjFMDyIkJAQSJKEs2fPeuyeGxckQAB444Dn7jkYSZLw/Vcy8JN3svDJ0XJvT4eIiIjILRiGByGEQHR0NHbu3Omxe8ZGBGLVZCPezSyB1Wb32H2dKa5rR0FVCwDgi9xKr86FiIiIyF0YhocQHR2NPXv2ePSe185LQFVzJ3bmu7/H8VCOlDQAAFL1wcgsqvfqXIiIiIjchWF4CLGxscjIyPDoPVdNNiIqWIN3Mr3bc/hYaSM0KgWumBmD0oZ2tJmtXp0PERERkTswDA8hLi4OJ0+e9Og9NSoFrpwdhy9zq1DfavbovXs7W9eGhMhATI4JAwCcrm712lyIiIiI3IVheAiRkZGwWCweD8TXzY+H2WbHh0dKPXrf3krq2xEfGYRkXTAARzgmIiIiGm8YhocghEBCQgI+++wzj953SkwYpsWG4Z3MEo/et7fShnbERQYiJlwLAChv7PDaXIiIiIjchWF4GLGxsdi1a5fH73vdvHgcK2tCblmTx+/d2mlFXasZ8ZGBiAhSI0ClQEVju8fnQURERORuDMPDSExMxMGDB0d8nc1mQ1ZWFjo6RreiumF2HDRKhVcepCvvCr5xEYEQQiAmXMuVYSIiIhqXGIaHER8fj7KyMlRVVbl8zY4dOxAXF4dly5YhJiYGW7duHfF9I4M1uHiqER8eKYPZ6tmewzUtjgf39CEBAIDocC0qGIaJiIhoHGIYHoZarUZSUhLeffddl84/evQo1q9fj0WLFuEnP/kJLr74Ylx//fUoKRl5/e/Vc+JR12rGzvzqEV87FnVdXSyigjUAAGOoFlXNnR6dAxEREZEnMAy7IDk5GVu2bBn2PJvNho0bN2LevHmYO3cuAGD69OlIT0/HAw88MOL7rphoQESQGh8eKRvxtWNR2xWGdV1hOCpYg/o277V5IyIiInIXhmEXTJgwAXv37oXdPnS5wuOPP466ujosXbq0z/srVqzAp59+ivLy8hHdV6NS4PIZMfgitxKtnZ7b9KK2xbEKHNkVhiODNGjusMLi5S2iiYiIiOTGMOwCk8kEhUIxZO1vfn4+/vrXv2L9+vVQqVR9joWFhSEtLQ1/+ctfRnzvDbPj0G6x4YvcyhFfO1p1rWaEB6qhVjq+HlHBagDg6jARERGNOwzDLhBCYNq0aXj++ecHPefmm2/GnDlzEBcX5/T47Nmz8c4774z43vOTIhEXEejRDThqW809JRLAuRXi+laLx+ZARERE5AkMwy6aMWMGtm7dipaWlgHHnn76aRQUFODCCy8c9PqUlBRUVlbi+PHjI7qvQiGwflYsduTX9JQvuFtdi7nn4Tng3IN0dV7cHpqIiIjIHRiGXaTX6xEXF4cnnniiz/tFRUV45JFHsGHDBqjV6kGvVyqVSE9PxyuvvDLie185JxY2u4QtR0dWczxata2d0IUMDMMskyAiIqLxhmF4BC688EI89dRTOHPmDACgvb0d69atw+zZs5GYmDjs9RMnTsSnn3464vtOjg7DJFMoNnmoq0RDmwURgb3CcBBXhomIiGh8YhgegdjYWCxYsADLli3D448/jpkzZ0KhUGDlypUuXZ+SkoITJ07AbB55qNwwJxaZRfUormsb8bUj1dxhRVjguYcAI4K6a4YZhomIiGh8YRgeoRUrVmDRokXYtGkTJk2ahGuuuQYKhWsfY0hICMLCwvDll1+O+L7rZ8YCAD7Odm+phMVmR7vFhlDtuZIPjUqBAJUCzR5s70ZERETkCQzDIySEwKxZs3DVVVdhwYIFLgfhbklJSaMqlUiICsKs+HB8muPeMNzc4Qi8Ydq+7eFCtWo0d7CbBBEREY0vDMMelpiYiN27d4/q2stmxCC7pNGtpRJN7Y7A23tlGHCE46YOrgwTERHR+MIw7GFJSUk4ceLEsLvZOXPZjBgAcGtXiZ6V4cC+YTg0UN1zjIiIiGi8YBj2sLCwMCiVShw9enTE1yZEBWFGXDi25FS4YWYOTR3dK8P9dtHTqlgmQUREROMOw7AXxMXFYdu2baO69rIZMcgqbkBJvXtKJboDb1i/MolQraqnhIKIiIhovGAY9oLo6Gjs3bt3VNdeNiMaAPDpUfesDje1O0oh+q8MhwawTIKIiIjGH4ZhL4iNjUVOTs6ork3SBWNabBg+cVPdcNMQK8MMw0RERDTeMAx7QWxsLAoLC0f1EB3gKJU4UtyA0oZ2mWd27gG6ECet1dotNlhso5szERERkS9iGPaCkJAQaLVaZGZmjur67q4Sn7phdbipw4KQABWUCtHn/e6yiRauDhMREdE4wjDsJdHR0aPuN5yiD8aUmDC3tFhr7rAO2HADONdqrYkdJYiIiGgcYRj2EoPBgEOHDo36+nXTo3G4uAFVzR0yzsqx6Ub/DTeAcyvDrBsmIiKi8YRh2EtMJhNyc3NHff2aqSZIErAtr0rGWQGtZiuCA5QD3g8JcITh1k6GYSIiIho/GIa9xGg0oqioaNTXT44ORUJUILYek7fFWmunDcEBA8skgjSOgNxmtsl6PyIiIiJvYhj2Ep1Oh8bGRtTW1o7qeiEE1kyJxu5TtbKu1raZrT3Bt7fugNxq5sowERERjR8Mw16iVCqh1+uxZ8+eUY+xZqoJZqsdO05WyzavNrMNwZohVoY7uTJMRERE44dLYVgIcakQ4oQQokAI8TMnx8OFEB8JIbKEEMeEELfJP9Xxx2Qy4eDBg6O+fkFyJCKC1NiaWynbnNrMNgQ5qRnuDshcGSYiIqLxZNgwLIRQAngawDoAUwHcIISY2u+0ewHkSpI0C8BFAP4shNDIPNdxx2AwICsra9TXq5QKrJpsxFfHq2TbDMNRJuFkZTiANcNEREQ0/riyMrwQQIEkSaclSTIDeBPAhn7nSABChRACQAiAOgBcQhyG0WjEiRMnxjTGJVNNaGy34GBh3ZjnY7NL6LDYndYMa5QKqBSC3SSIiIhoXHElDMcBKO71uqTrvd6eAjAFQBmAowB+KEnSgKVKIcQdQogMIURGdbV8da7nK5PJhJKSkjGNsWKiAQEqBbYeG3upRFtXCYSzmmEhBAI1Sq4MExER0bjiShgWTt6T+r1eC+AIgFgAswE8JYQIG3CRJD0nSdJ8SZLmGwyGEU51/AkNDYXNZhtTi7UgjQrL0vT4IrcSktT//5aRae8KuoFOVoYBR0huY80wERERjSOuhOESAAm9XsfDsQLc220A3pccCgCcATBZnimOX0IImEwm7N27d0zjrJlqQmlDO3LLm8Y0TmtXGHa26QbgqBtu5cowERERjSOuhOGDANKFECldD8VtBLC53zlnAawGACGECcAkAKflnOh4ZTQax7QtMwCsnmKCEMCXuWPbja571dfZA3RA18owa4aJiIhoHBk2DEuSZAVwH4DPAeQBeFuSpGNCiLuEEHd1nfZrABcIIY4C2AbgYUmSatw16fFEr9cjJydnTGMYQgMwMz4C20+MNQw7Vn2dPUDX/T5XhomIiGg8cb4E2I8kSVsAbOn33r96/boMwCXyTs0/GI1G7N69e8zjrJxkwJPb8lHb0gldSMCoxjgXhgdZGQ5Qoaq5Y9RzJCIiIvI13IHOy4xGI0pKSsb88NuqyUZIErAjf/RdOrpLIIZaGeYOdERERDSeMAx7WXBwMJRKJc6cOTOmcabHhkMfEoCvjo8+DPc8QDdUzTDLJIiIiGgcYRj2AXJ0lFAoBC6aZMCOk9WwjnI3uvbuB+gG6SYRqFFyO2YiIiIaVxiGfYAcHSUAYOUkIxrbLThc3DCq61uHeYAuOMCx6cZYSzqIiIiIfAXDsA+Qo6MEACyfqIdSIbD9+Oi6SrSZbRAC0KoGqxlWwWaX0Gkd3cozERERka9hGPYBRqMR+fn5Yx4nTKvG/KRIfDXaMNxpRZBaCYXC2aaD51aM21k3TEREROMEw7APMBgMKCsrg90+9hXXVZONOF7RjPLG9hFf22q2IXCQh+cAQKt2hOEOK8MwERERjQ8Mwz4gMDAQAQEBOHXq1JjHWjnZCAD4+sTIu0q0m62DbsUMAIFqrgwTERHR+MIw7CNMJhP27ds35nHSjSGIiwgcValEm9nWE3id0aodX5cOC2uGiYiIaHxgGPYRJpNJlo4SQgismmzE7oIadI6wnKHdYusphXCm+1i7hSvDREREND4wDPsInU6Ho0ePyjLWyskGtJltOHimfkTXdVrsPau/znSH4U6GYSIiIhonGIZ9hMFgQEFBgSxjLU7VQa0U2DnCrZk7rEOvDAdyZZiIiIjGGYZhH9HdUcJmG3vQDNKoMD8pCjvya0Z0XYfFNmiPYaBXNwnWDBMREdE4wTDsIwICAhASEoKTJ0/KMt7yiXrklTehqrnD5Ws6LHYEDrL7HMCVYSIiIhp/GIZ9iFwdJQBgRboBALC7wPXV4Q6LbZiaYUXPeURERETjAcOwDzGZTMjMzJRlrKkxYYgK1mDnyZGF4YChyiQ0yp7ziIiIiMYDhmEfEhUVhWPHjskylkIhsCxNjx35NZAkyaVrOqz2oVurqRiGiYiIaHxhGPYhRqMR+fn5so23YqIBNS2dOF7RPOy5NrsEs3Xo1mpqpYBSIVgzTEREROMGw7AP0ev1qKyshNlslmW85el6AHCpxVr3Bh1DrQwLIaBVKdhNgoiIiMYNhmEfolarER4ejtzcXFnGM4VpMckUih0u1A13B1ytauivRKBGyZVhIiIiGjcYhn1MdHQ09u7dK9t4y9P1OFBYh3bz0AG2uw54qJVhAAhQKVkzTEREROMGw7CPMRqNOHTokGzjLZ9ogNlqx4HCuiHPczUMB2oYhomIiGj8YBj2MTqdDjk5ObKNtzA5ChqVAjtPDl033FMmMcQDdN3HWTNMRERE4wXDsI8xGo04ffq0bOMFapRYmByFncNszdzR9QBdwHArw2rlsCUXREREROcLlbcnQH3pdDrU1dWhvb0dgYGBsoy5LF2P3396HFVNHTCGaZ2e01MmMcSmG4CjjKKl0zrsPSVJwpd5VXjzwFlUNHVgcnQYbl+RgsnRYSP/AYiIiIjchCvDPkapVEKn08m2LTMALJ3gaLG293TtoOd0h+FAzfBheLiVYZtdws/eO4rbX83A8YpmGEIDsPVYBa74+y68k1E8wtkTERERuQ/DsA8yGo2ydpSYGhuG8EA19hQMFYZdrRlWotM6dM3wrz46hrcyinHPRRPw9U8vwr9vW4gdD63E4lQdfvpuNj7LqRj5D0FERETkBgzDPshkMmH//v2yjadUCCxOjcKe04PXDbtaJhGoVgy5MvxJdjle2VuEHyxLwUOXToZa6fiKRQZr8MKt8zE7IQI/eScLZQ3to/hJiIiIiOTFMOyDEhISkJWVJeuYF0zQo7iuHcV1bU6Pn1sZHv4Buu6H7fpr7rDg8Y+OYWZ8OB5eN3nAca1aiX/cMAc2u4RfbpKvYwYRERHRaDEM+6DY2FiUlZWhpaVFtjEvmKADAOw55Xx1+Fyf4eHLJAZbGX5qewFqWjrx6w3Te1aE+0uICsKDF6dj2/Eq7D01eNkGERERkScwDPsgtVoNo9GIL7/8UrYx04wh0IcEYM8gAbR7tXfYHei6aobtdqnP+/WtZry2twgbZsViVkLEkGPcekEyTGEB+L+tJyBJ0pDnEhEREbkTw7CPSkhIkDUMCyFwwQQd9pyqdRpAu8skAlRDfyW6j5ttfR+ie3lPIdrMNtyzMm3YuWjVSty7Mg2ZRfXILKp39UcgIiIikh3DsI+KjY0dVXu1qqoqfPDBB2hoaBhw7IIJOlQ3d6KgamD5RafFhgCVAkKIIcfvDsO9O0p0Wm34z74iXDzFhImmUJfmee28eIRqVfj3nkKXziciIiJyB4ZhH5WcnIzc3FxYLBaXr3n66aeRkpKCe++9F8nJydi+fXuf40vTHP2GnZVKdFhsw5ZIAOd2qOvs9RDdF7mVqGs146bFiS7PNUijwsYFCfg0pwIVjR0uX0dEREQkJ4ZhHxUWFobw8HB89tlnLp3/1ltv4eGHH8aNN96IO++8E6tWrcJ1112Hurq6nnMSooIQHxno9CG6dott2IfngF5lEr1Wht86WIy4iEAsTze4NNdu31mUBJtdwqYjpSO6joiIiEguDMM+LC0tDe+///6w55WVleGOO+7AVVddhfj4eADArFmzEB0djUceeaTPuRdM0GHf6TrY+j0A12Gxu7Yy3K9MorShHTvza3D9/AQoFUOXWPSXrA/GnMQIbDrMMExERETewTDsw1JTU/H1118Pe96tt96KSZMmIS2t78NrK1aswOuvv96nRdsFE/RobLcgr7ypz7kdFhsCRxKGux64+/RoOQDgyjmxw17rzFVz4nC8onnAfIiIiIg8gWHYhyUmJqKqqgrHjh0b9Jy33noLBw4cwOrVqwccM5lMMBqNePbZZ3veW9LVb3h3Qd9SiQ6rvaceeCgBqr41w1uOlmNabBiSdMHD/0BOXD4jBiqFwOasslFdT0RERDQWDMM+TKVSYfr06XjyySedHm9tbcX999+PdevWQavVOj1n5syZeO2113pem8K0mGAIxu5+D9F1WGzQDtNWDehbM1zW0I5DZxtw2YwYV3+kAXQhAVicqsMXuZWjHoOIiIhotBiGfdyMGTPwwQcfwGq1Djh29913w2g0YsqUKYNeP3nyZOTl5aG29lz4XTJBh8zCOlhtvdujubgyrD5XM/xZTgUAYN30aJd/HmcunmJEQVULztS0jmkcIiIiopFiGPZx8fHxCAoKwt/+9rc+72/duhXvvfce1q5dO+T1Wq0W8fHxeOONN3reW5SiQ6vZhmNl5+p0u/sMD0ej7C6TsGP7iSqkGUOQaggZwU800OopJgDAtjyuDhMREZFnMQz7OCEELrroIjzxxBMoL3c8rHbq1Cls3LgRl19+OcLCwoYdIy0tDZs3b+55vSglCgCw/8y51WKzzQ6NK2USXSvD9W1m7D9dh4smjqydmjMJUUGYHB2KrSyVICIiIg9jGD4PpKamYvr06Zg3bx5uv/12zJ8/H/Pnz8e0adNcuj4tLQ0ZGRk92zAbw7RI0QfjwJlzPYjNVjsClK7XDO84WQ2zzY4LJ409DAPAxVNMyCisQ1OH65uMEBEREY0Vw/B5YvXq1Vi2bBny8vKwYcMGLF261OVrDQYDOjs7+3SlWJQShQNnzvUbdtQMuxKGHWUSW3MrEahWYkFy1Ah/EueWpulhl4D9p+uGP5mIiIhIJgzD5wkhBKZNm4Y1a9YgJSVlxNcmJSXh448/7nlvYUoUmjqsOFHRDMCxMqwZwcqw2WrHkgk6lzbqcMXcpAho1YoBLd+IiIiI3Ilh2E8kJCTgm2++6Xm9KNXRb7i7bthsda1muPc5y9P1ss0vQOVYZWYYJiIiIk9iGPYTiYmJyMrK6nkdFxGI+MjAnrKETqutpwRiKL07TizuCtRyWZqmR35VC6qaOmQdl4iIiGgwDMN+Ijo6GjU1NaipObfyujAlCgcK62Cx2WGX4NLKsKpXKcUkU6isc1yW5lhp3n2Kq8NERETkGQzDfkKlUsFkMuGrr77qeW9xig51rWbkdvUbdiUM96ZQCFnnODUmDKFaFQ6cqZd1XCIiIqLBMAz7kZiYGOzevbvn9cKufsM786sBwKVNN6qa3VfCoFAIzE2MRGYRO0oQERGRZzAM+5Ho6GgcOnSo53WSLgimsADszHeUJbiyMnywa9VW7hKJbvOTInGysgWNbew3TERERO7HMOxHYmNjkZ+f3/NaCIFFKTrs79p8w5XWage6uk+kmca2BfNg5iVFAgAOFbNUgoiIiNyPYdiPGAwG1NXVoa7uXBlCd6kEAAS40DM4o8gRUrt3s5PbrIQIKBUCmYUMw0REROR+DMN+RKVSwWg04uuvv+55b3HquTA83Mpwh8WG412bdHRa7G6ZY3CAClNiQpHBumEiIiLyAIZhPxMdHY0DBw70vJ5gOFfuMNwDdMfKGvts3+wu85OikFXcCIvNffcgIiIiAhiG/Y7BYEB2dnbPayEEdMEaAMOH4SPFjQCACYZgdFptbpvjnMQItFtsOFnZ7LZ7EBEREQEMw37HaDT2eYgOAGbGhwMA6trMQ16bVdyAmHAtEqKC3LoyPDM+AgCQU9rotnsQERERAQzDfsdoNKK0tBR2+7kwO6MrfB4dJnxmlzRgZnw4AlQKmN0YhpOighAaoEJ2CcMwERERuRfDsJ8JDg6GEAInTpzoeW+CIRgAcHSI8NnQZkZhbRtmJUQgQKV068qwQiEwPS582HBORERENFYMw35GCIHo6Gjs3bu35z17V5u0ocJwVtex2fER0KgU6LS4r2YYcJRuHC9vdusKNBERERHDsB8yGAx9dqLrDpzNnVa0dlqdXpNd3AAhgOldZRLuXBkGgOlx4TDb7HyIjoiIiNyKYdgPGQwG5OTk9Lzuvfp6pLjB6TVHSxuRogtGmFaNAJXS7Su23Q/1sW6YiIiI3Ilh2A8ZjUacPn2653XvVd6MQXZ+y6towpSYMABAgNr9K8OJUUEI06pYN0xERERuxTDsh4xGIyoqKmC1OkoiuoNtqj7Y6c5vzR0WFNe1Y0pMKADHTnVmmx12u3u2ZAYctc0z4sNxtLTBbfcgIiIiYhj2Q1qtFkFBQcjKygJwrkxiyQQdDp9t6Nllrlv3Fsy9V4YBwOzmHeKmxoThZGULrNyJjoiIiNzEpTAshLhUCHFCCFEghPjZIOdcJIQ4IoQ4JoT4Rt5pktxMJlPPtsydVjs0SgUWJEehpdOK4xVNfc7NK3e87gnDKqXjOot7Q+qk6DCYrXYU1ra59T5ERETkv4YNw0IIJYCnAawDMBXADUKIqf3OiQDwDIBvSZI0DcB18k+V5KTX63HkyBEAjpVhjUqBeUmRAIDMor51w3nlTQgPVCMmXAsA0CiF4zo3r9hOjnaUZfQP50RERERycWVleCGAAkmSTkuSZAbwJoAN/c65EcD7kiSdBQBJkqrknSbJzWAwIDc3FwBgttkQoFIgPjIQprCAAQ/R5ZY3Y0pMKIRwhGCNyjNlEmnGECgVAsfL2V6NiIiI3MOVMBwHoLjX65Ku93qbCCBSCPG1ECJTCHGLXBMk9zAYDCgsLATgKHfQqBQQQmB+chQyCs89RGezSzjRq5MEAKiVjq+Nxc0dJbRqJVL1wVwZJiIiIrdxJQwLJ+/1byOgAjAPwOUA1gL4pRBi4oCBhLhDCJEhhMiorq4e8WRJPgaDAZWVlbBarTDb7D2rvfOTIlHW2IGyhnYAQGFtKzos9j5huPtciwcebJsUHdrzAB8RERGR3FwJwyUAEnq9jgdQ5uSczyRJapUkqQbADgCz+g8kSdJzkiTNlyRpvsFgGO2cSQYBAQEIDg5GdnY2zFY7AnrCcBQAIKOrbrj74bmpTlaG3d1rGHA8tFdS346mDovb70VERET+x5UwfBBAuhAiRQihAbARwOZ+53wIYLkQQiWECAKwCECevFMluZlMJuzfv9/RTaIrDE+JCUWQRonMrlKJ4+XNUCoE0owhPdd5cmW4+yG6k1wdJiIiIjcYNgxLkmQFcB+Az+EIuG9LknRMCHGXEOKurnPyAHwGIBvAAQAvSJKUM9iY5BsMBgOysrIc3SS6VntVSgVmJ0TgYNdDdPlVzUjSBUGrVvZc132uu7dkBoDJXSvSLJUgIiIid1C5cpIkSVsAbOn33r/6vf4TgD/JNzVyN71ej9zcXETPtvf0DgYcdcNPbS9Aa6cV+VUtSO+1Kgz0eoDO5r4d6LrFhmsRqlXxIToiIiJyC+5A58cMBgNOnz6NTqutp/QBAOYkRcIuOfoNF9W2Id0Y2uc6T5ZJCCGQbgxBfmWL2+9FRERE/odh2I91d5Ro7zT3CcOz4yMAAB8cLoXNLvWpFwYAddemG554gA5w9Bs+Vc0wTERERPJjGPZjGo0GoaGhaKwo7ukmAQCRwRqk6oPxweFSABgQhgM8uDLcff+aFjMa2sweuR8RERH5D4ZhP2cymVBXcqrPyjAAzE6MAAAIAUwwOK8Z9sQDdMC5MF5QxdVhIiIikhfDsJ8zGAxoqTjTZ2UYAOYkRgIAJAkI1Cj7HPNkzTAApBkcNcsMw0RERCQ3hmE/p9fr0VZV1KebBADMSYgY9Jpz3SQ8E4bjIgOhUSlYN0xERESyYxj2c5cltCJ77VE8dmgp8NfpQPbbAAbWCffmyR3oAECpEEjVB3NlmIiIiGTnUp9hGp+mS3lYH5EBjegKtY3FwEcPAADKYy8f9LpzD9C5v89wtzRjCLJKGjx2PyIiIvIPXBn2Y6uxCxph7fumpR3Y9ivkV57b8a3DYutziqcfoAMcYbikvn3AXIiIiIjGgmHYj4VjkC2OG0uQ36sk4VhZ393flAoBpUJ4rGYYcIRhSQLrhomIiEhWDMN+rBGhzg+Ex+N0dSuEY28NHD5bP+AUtdLzYRhgRwkiIiKSF8OwH9uGZTD3LxtXBwKrH0VhbSsWJkchLiIQh4sbBlyrVio89gAdAKTog6EQwCmGYSIiIpIRw7AfyxFT8BHWoNYaBLsElNj1KFn2B2Dm9SisaUWKPhhzEiNwuGjgynCASuHRleEAlRJxkYE4U9vmsXsSERHR+Mcw7OdyxBT8Fd+H5rdtWCs9iccLp6Kpw4LaVjOS9cGYkxiJssYOVDR29LlOrVR49AE6AEjWBaOwptWj9yQiIqLxjWGYoFarEREejnWxVnyZV4VNh0sBOMLnnK5tmY8U910d1nh4ZRhwlEoU1rRCkjzX0o2IiIjGN4ZhAgCYTCYY285AF6zBox8eA+AIn9Niw6BRKnD4bEOf89VKhUf7DANAki4YzZ1W1LWaPXpfIiIiGr8YhgmAY1vm3KNZuPuiCT3vJemCEKBSYmps2IAwrPHwA3QAkKIPAgAU1rJUgoiIiOTBMEwAAIPBgGPHjuE7i5J63tOqlQCAOYkRyC5t6FMWofZCmUSyLhgAcKaGD9ERERGRPBiGCQAQHR2NU6dOIVCj7Hmvu7/w3MRIdFjsOF5+bpMOjVJ4/AG6+MggKBWCD9ERERGRbBiGCQCg0+nQ1taG4uJiaLq2W/7HVwUAgNkJEQCAIyUNPed74wE6jUqBuIhAlkkQERGRbBiGCQCgUCgQFxeHD7d8DrPNjqhgDb46XoWjJY2IjwxEVLAG2b0233A8QOfZMAwAyfpghmEiIiKSDcMw9YiNjcW2HXsAAI9eMRVhWhX+8VU+hBCYFR+OrN4rw154gA4AUnRBKKxpY3s1IiIikgXDMPWIiYlBdlYWAGB6XDi+tywFW3MrkVfehJnxEcivakFLpxWAdx6gAxzt1Vo6rahpYXs1IiIiGjuGYeqxPrkDX60pwOmAGzHhv4twR0QmQgNUeOqrAsxOiIAkATmljQCAAKUCZi+E4RS9o6NEEUsliIiISAYMwwQAmC7l4dtBu5EUDigEIBpLEPT5j/C79OPYklOO4AAVACCrq25YrVTAYvV8qUKyvru9GsMwERERjR3DMAEAVmMXNLD2fdPSjnVVz0GjVOC9zBLERwYiu8SxMqxWCa+UScRHBjraq3FlmIiIiGTAMEwAgHA0O31f0VSKb82KxQeHSxETrsWRrpVhjVLZp8/wnj178NRTT6G8vNyt81QrHe3Vzta1u/U+RERE5B8YhgkA0IhQp++X2nV4J7MEZpsdBwvrUdrQjpqWTqhVAmabHXa7HbfddhsuueQSPPnkk5gyZQr27Nnj1rkmRAWiuI670BEREdHYMQwTAGAblsEMVZ/3JHUg2pb/DxalRPV5v3znq7j3yJXIU25E0/8mIuDEh7jrrrtw0003YcWKFbjmmmvQ3u6+lduEyCCU1DMMExER0dgxDBMAIEdMwUdYgwaEwi4BDVIoxPq/Y+Ka7+PNOxbjp2snAQC+pdiFCfv+B2GdFVAIIEI048mLgSUhJQCAuXPnIiQkBI899pjb5poQFYSaFjNaO63Dn0xEREQ0BIZh6pEjpuBJcTuuP7YGF32oB2ZeDwAQQuDelWlYkqrDQ6q3EST69vgNUNiwGrt6zl22bBn+/e9/w2azuWWe8ZGBAICSetYNExER0dgwDNMA6enpOH78OKqrq/u8f/dFExAnapxe0/sBvKSkJCiVSrz++utumV9iVBAAsG6YiIiIxoxhmAbQarVITU3FP//5zz7vL+/YjsE6C/d+AE8IgRkzZuCVV15xy/wSusMw64aJiIhojBiGyak5c+bgxRdfhCSdi79i26+gEAPPleB4AK+3qVOnYu/evTCb5d82WResQaBaiWK2VyMiIqIxYhgmp9LS0tDa2orXXnvt3JuNJYOenyOm9HkdGRmJ8PBwfPzxx7LPTQjhaK82hpXhwppWfJLt3p7IRERE5PsYhskphUKBSy65BD/5yU96aofb1FFOzx2sR3FKSgo2bdrklvklRAaNqWb4wbeO4N7XD+FUdYuMsyIiIqLzDcMwDWrSpEmYMGECZs6cibVr1+Le96vQaVf2OccM1YASiW6pqalu24AjIcoRhnuXcYxE9056R842yDcpIiIiOu8wDNOQ1q5di9WrV0OhUCB8+Q/wseISNCAUEoAGhOIjrBlQItEtMTERxcXFqKlx3oFiLBKigtBqtqG+zTLia5s6zl1zlh0piIiI/Jpq+FPInwkhkJ6ejvT0dABADiKRA+fhtz+NRoOYmBh8+umnuPnmm2WdV0JXr+HiujZEBWtGdG1pr/7ElU0dss6LiIiIzi9cGSa3iouLw1dffSX7uGNpr1be2N7r1wzDRERE/oxhmNwqPj4emZmZso/bE4ZH0V6toau0IkUfjLpW+Vu/ERER0fmDYZjcKj4+HqdOnYLdbpd13JAAFaKCNaOq+W1sd4ThhKignl8TERGRf2IYJrcKCwuDRqPBwYMHZR87PjIQJaMok+gOwPGRgQzDREREfo5hmNwuLi4OX3/9tfzjRgSOqua3sd2CkAAVdMEaNHVYYLePrj0bERERnf8YhsntjEYjDhw4IPu4sRGBKK1vH3Gv4aZ2K8K0KoQHqiFJQHOnVfa5ERER0fmBYZjcLiYmBnl5ebKPGxsRiHaLreeBOFc1tlsQFqhGeKAaANDEUgkiIiK/xTBMbhcdHY2zZ8/K/hBdXIQWAFDaMLKOEk3tFoQHqhHWFYZZN0xEROS/GIbJ7UJDQyFJkuyrw7ERjo03ykYahjscK8NBGsfW0u0Wm6zzIiIiovMHwzC5nRACsbGx2Llzp6zjxo0yDLearQgJUCFI49iAsZU1w0RERH6LYZg8wmQyyd5eLSpYgwCVAmUj7CjRbrZBq1YiOEDZ85qIiIj8E8MweYTJZMLRo0dlHVMIgbiIwBHXDLeZbQjSKBGk7loZZhgmIiLyWwzD5BEmkwmFhYWyjxsbETiiMglJktBu6QrDXSvDbWaWSRAREfkrhmHyCJ1Oh/r6erS0tMg6bmyEFqX1rofhTqsdkgRHmURXzXAbV4aJiIj8FsMweYRKpUJkZCT2798v67ixEYGoau5Ep9W1QNtdHxykUUKrVkAIoI0P0BEREfkthmHyGKPRiIyMDFnH7G6vVtnY6dL5bZZzYVgIgSC1kjXDREREfoxhmDxGp9MhKytL1jHju8Kwqw/RtXfVB2vVjnrhoAAVyySIiIj8GMMweYzRaMSJEydkHXOkG2+0mx274HX3GA7SKPkAHRERkR9jGCaPMRgMKC4ulnXM6HDHlsyuhuHu4BvYvTKsUaG1kyvDRERE/ophmDymu6NEc3OzbGNq1UroQwJQ1ujiynBXzXBg11bMWrXC5YfviIiIaPxhGCaPUalUiIqKkr2jRFyEFqUNru1C17ubBAAEqBTotNhlnQ8RERGdPxiGyaMMBoNbOkqU1re5dG73w3LdZRJatZIrw0RERH6MYZg8Sq/Xy74tc2xEIMobOyBJ0rDntlsGrgx3cGWYiIjIbzEMk0cZDAbZO0pEh2nRZrah2YXNM7rLJAJ7wjBXhomIiPwZwzB5lNFoxNmzZ2Uds7ujREXj8HXDPQ/Qqc89QMeVYSIiIv/FMEweFRUVhfr6erS2tso25kjCcKfVBqVCQKV0fPW5MkxEROTfGIbJo1QqFSIjI3HgwAHZxowOG0EYttgRoDr3tefKMBERkX9zKQwLIS4VQpwQQhQIIX42xHkLhBA2IcS18k2RxhuTySRrRwlTVxgudyEMm212aHqF4e6VYVceviMiIqLxZ9gwLIRQAngawDoAUwHcIISYOsh5fwDwudyTpPFFr9cjKytLtvE0KgX0IRpUNLkQhq12aJR9V4btEmC1MwwTERH5I1dWhhcCKJAk6bQkSWYAbwLY4OS8+wG8B6BKxvnROOSOjhKmMC0qXNiFrtNqR4C678owAHRYWDdMRETkj1wJw3EAinu9Lul6r4cQIg7AVQD+NdRAQog7hBAZQoiM6urqkc6VxgmDwYCioiJZx4wJ16KiqXPY8/qvDHcH404r64aJiIj8kSthWDh5r//fKf8NwMOSJA25vCZJ0nOSJM2XJGm+wWBwcYo03nR3lGhrc23XOFeMZGVY07UaDABargwTERH5NVfCcAmAhF6v4wGU9TtnPoA3hRCFAK4F8IwQ4ko5Jkjjj0qlQlRUFA4ePCjbmDHhWtS3WYYNtZ1WW59uElwZJiIi8m+uhOGDANKFEClCCA2AjQA29z5BkqQUSZKSJUlKBvAugHskSdok92Rp/DCZTLKG4e6OEpXDPERntg7sJgFwZZiIiMhfDRuGJUmyArgPji4ReQDeliTpmBDiLiHEXe6eII1POp0O2dnZso0XEx4IYPhew2abnSvDRERE1EPlykmSJG0BsKXfe04flpMk6btjnxaNdwaDAcePH5dtvOjwAAAYtr1ap8UOXXCv1mpdK8Od3HiDiIjIL3EHOvIKg8GAwsJC2caLHsHKsMbJynAHt2QmIiLySwzD5BVyd5QICVAhJEA17C50A1qrdQXjTtYMExER+SWGYfKK7o4Scm7LHB2uHXZl2NFN4lxrte4wbLZxBzoiIiJ/xDBMXmMymXDgwAHZxosO0w5bM9y/m4RG6QjGFj5AR0RE5JcYhslrdDodsrKyZBvPlZXh/mFYrXLsKWO2MQwTERH5I4Zh8hqDwYC8vDzZxosO06K6pRPWIYJtp7VvazV1V/2whWGYiIjILzEMk9cYDAYUFRXJNl50uBY2u4SaFrPT43a7BKtd6lsm0V0zzDIJIiIiv8QwTF6j0+nQ0NCA5uZmWcaL7tqFrryx3enx7lKIvjXDij7HiIiIyL8wDJPXKJVK6HQ67NmzR5bxosO7t2TudHq8e2ON3t0kesokrOwmQURE5I8YhsmrYmNjsXv3blnGMoV1h2HnD9F12hy9hHuvDCsVAkqFYM0wERGRn2IYJq8ymUyy9RrWBWugVopB26t11wUHKPt+7dVKwTIJIiIiP8UwTF4VExMjW0cJhULAGKpF5SDt1Tq7w7C6fxhW8AE6IiIiP8UwTF4VHR2N0tJSmM3OO0CMlCksYNiVYU2/leEAlYJlEkRERH6KYZi8KiAgAOHh4di3b58s40WHawetGe4JwyquDBMREZEDwzB5XWxsLHbs2CHLWMZQ7eDdJKwDu0kAjjDMlWEiIiL/xDBMXhcdHY39+/fLM1a4Fi2dVrR0WgccG2xlWKNSwGJjazUiIiJ/xDBMXhcXF4fs7GxZxureeKPCyUN03au/aqXo875aqehZNSYiIiL/wjBMXhcXF4eKigrU19ePeayheg2be8Jwv5VhJfsMExER+SuGYfI6lUqF2NhYbNmyZcxjnduFbmAYtnaVQgwIw+wmQURE5LcYhsknJCQk4IsvvhjRNXa7HS+88AIeeeQRFBUVAXC0VgPgtL3aUGUS7CZBRETkn1TengAR4AjDe/fudfn8trY2rFixAiUlJTAajfjnP/+JL7/8EvPnz0eoVuV04w3LIGUSaqUCrU4euCMiIqLxjyvD5BOSkpJQWFiIurq6Yc+VJAlXXHEFWltbcfvtt+Oaa67B4sWLcd1118FqtSI6TOt0ZdhqH7xMwsxuEkRERH6JYZh8QkBAABITE/Hmm28Oe+6TTz6Jo0eP4uqrr4ZK5fjLjcWLF8NqteJvf/sbosO1qHDSa7h7ZVjVr0xCo1TAbLXJ8FO4V0unFZlF9fgitxJf5FYip7QRHRbfnzcREZEvY5kE+YwJEyZg06ZNuOeeewY9p7S0FI899hiuueYaaDSanveFEFi2bBmefvppXPPHNSioqhlwbXcvYbWif5mE8Nk+w2arHZuOlOK9zBIcKKyD1G+agWollqbpcdPiRFw40QAhhPOBiIiIyCmGYfIZkyZNwosvvoiOjg5otVqn59xyyy2YPHkykpKSBhxLT0/HJ598grbS46hqDoXNLkGpOBcOe2qGVf1Whn20m8TWYxX41ce5KKlvR6ohGPetTMOs+AiYwrSwSxJK6ttx4EwttuRU4MuXKzErPhy/2jAdsxIivD11IiKi8wbLJMhnREVFwWAw4MUXX3R6/J133kFGRgZWrVrl9LhSqcSsWbOQ9eUHsNkl1Lb0LZWwdpdJDFgZ9q1uEm1mK3701hHc8VomQgJUePm7C7Dt/12IH18yCRdPNWFGfDhmJUTg8pkx+N8N07H74VX447UzUdbYgaue2Y2ntxfAbvfNlW4iIiJfwzBMPmXWrFl47rnnBrzf2tqK++67D+vWrRt01RgAJk6ciOMZuwAMbK/WUybhrLWaj6wMVzR24Ppn9+LDI6X44ep0bL5vGVZONg5Z/qBRKXD9/AR89eMLccXMWPzp8xO48z+ZrCcmIiJyAcMw+ZTp06ejsLAQH3/8cZ/3b775Zuj1ekyZMmXI6+Pj49HR0ghLQ8WALZktNjtUCjEgWAb4SJlEWUM7rv3XHhTWtOGFW+fjR2smQqNy/V/RUK0aT26cjcfWT8WXeZX47ssH0MKWcURERENiGCafolarsWrVKtx555092zM/+uij2L59O6644ophr1coFEibOAkdpzNQ2dyvTMIuDegkAfhGmURlUwdufH4fGtsteP32RVg12TSqcYQQuG1pCv56/WwcLKzHHa9moPM86JRBRETkLXyAjnzOrFmzUFZWhokTJ8JkMqGsrAw333wzgoKCXLo+JSkBp/KODth4w2y1D+gxDDjCsF3CgAfuPKXNbMVtLx9EdXMnXvvBIsyMjxjzmFfOiYNdkvD/3s7CT97JxpPfng2FF342IiIiX8cwTD5HCIF169bh7NmzaG1txYYNG/q0URvO1WlWPJl8FIl7FwG58cDqR4GZ18Nqdx6Gu0sRzFY7AjVK2X4OV0iShJ++k428iia8dOsCzE2MlG3sq+fGo7KpE3/47DimxITinovSZBubiIhovGAYJp8khHDaPm0406U8rA/eA0d0loDGYuCjBwAAFuvkAQ/PAeceqLPY7QiEZ8Pw8ztP45Oj5fjZuslYOdko+/h3XZiK3PIm/N/nJzAnIRJLJuhkvwcREdH5jDXDNK6sxi5o0O+hMUs7sO1XsNjtA9qqAYCqq3zA6uGNN3LLmvCnz0/gkqkm3Lki1S33EELgd1fPQIo+GPe/cRh1rWa33IeIiOh8xTBM40o4mp0faCyB1SY57c6g6iqdsHqwo0SHxYYH3zqMiCANfn/NTLfuHBcSoMLT35mLxnYzHv0wx233ISIiOh8xDNO40ohQ5wfC43taq/V3rkzCcyvDf/syHycrW/Cna2ciKtj1eujRmhwdhh+uTsfH2eX49Gi52+9HRER0vmAYpnFlG5bB3L8UXh0IrH4UFpvUswrcW3fphKdWhvMrm/HCztO4dl48Lpokf53wYO66cAJmxIXjlx/moKnD4rH7EhER+TKGYRpXcsQUfIQ1aEAo7BJQbAlH85q/ADOvh8Vmh8bJA3TdvYctHqgZliQJv9iUg+AAFR5ZN9nt9+tNpVTgd1fPQG2rGX//Mt+j9yYiIvJVDMM07uSIKXhS3I45W6Zh1v5Lscl2AQDAarcPvTJsd//K8IdHyrD/TB0evnQydCEBbr9ff9PjwrFxQSL+vacQBVWD1FcTERH5EYZhGreSYw1QVOdjc1YZAMfKr7PWat0rw+7uJtFhseGPnx3H9LgwbFyQ4NZ7DeUnl0xEkEaJ//0o12tzICIi8hUMwzRuxcbGwlJ2AgcL61Ha0A6LbbAd6LrCsJsfoPvPviKUNXbgkXVTvLobnC4kAD+8eCJ25tdgd0GN1+ZBRETkCxiGadzS6/Voa6iGZDXjo6wyWG2S024SnniArqnDgqe3F2B5uh5L0/Ruu4+rblqciNhwLf70+QlIkmf7KxMREfkShmEat1QqFXQ6HeKlamw+UjboyrAnHqB7fsdp1LdZ8PClnn1objABKiUeWJ2OI8UN+DKvytvTISIi8hqGYRrXTCYTjG2FyC1vwvGK5kHKJNz7AF1jmwUv7TqDy2fEYHpcuFvuMRrXzItHsi4If956AnYP9lgmIiLyJQzDNK4ZDAZ0luahe4M3pw/QuXk75lf2FqLVbMN9q9LcMv5oqZUK/PDidByvaMa241wdJiIi/8QwTOOayWRC0amTWJKqAwCnrdW6V4YtbqgZbjNb8fLuM1g92YgpMWGyjz9W62fGIi4iEP/65pS3p0JEROQVDMM0rhmNRhQXF+PyGTEAgFPVLQPOUbmxm8QbB4pR32bBPSt9a1W4m0qpwO3LU5BZVI+DhXXeng4REZHHMQzTuBYeHg6LxYIpoWYAwOGzDQPO6e4mIffKsNlqx/M7TmNxahTmJUXKOracvr0gEVHBGvzra64OExGR/2EYpnFNCAGTyYS8IwcGPUftpk03Ps0pR0VTB+5cMUHWceUWqFHi1iXJ2Ha8CicruSsdERH5F4ZhGveMRiMyMzN7Xuf3C3wqN3WTeGVPIVL0wbhwokHWcd3hliVJCFAp8OreQm9PhYiIyKMYhmncMxgMyMnJ6Xn9ydHyPsfVCvn7DB8tacShsw24eXGSV3ebc1VksAbrZ8Xi/UOlaOqweHs6REREHsMwTOOe0WhEwalz9bCfHq3oc1zZ01pNvpXhf+8pRJBGiWvnx8s2prvduiQZbWYb3s8s8fZUiIiIPIZhmMY9g8GA8rIySJIEY2gATlQ2o6DqXFeJc2US8qwM17Z04qPsMlwzNx5hWrUsY3rCjPhwzE6IwKv7irhFMxER+Q2GYRr3goODoVSpYGuuxfpZsQCAz3LOlUqoZW6t9lZGMcxWO269IEmW8Tzp1guScLq6FbsKarw9FSIiIo9gGCa/YDRGw1JThPjIQMxPisQnvUoluluryVEmIUkS3j5YjIUpUUgzho55PE+7bEYMooI1ePNAsbenQkRE5BEMw+QXdEYjLDVFUCkVWDcjBnnlTThT0wrg3MqwHA/QHSysR2FtG66fnzDmsbwhQKXElbPj8EVuJepbzd6eDhERkdsxDJNfiNQZYK45C7VCYN30aACOPsCAoxexUiFkaa321sFihASocNmM6DGP5S3XzY+H2WbHh0dKvT0VIiIit2MYJr+wPrkDx9Zk4NtbZiL25QV4wHAYW3q1WFMpxJg33WjusGDL0XKsnxWDII1qrFP2mikxYZgRF463M9hVgoiIxj+GYRr3pkt5uEN3EEmhdghIQGMx7m/9B1LLt6C4rg0AoFYqxlwm8Ul2OdottvO2RKK36+bHI7e8CTmljd6eChERkVsxDNO4txq7ECCsfd5T2zvwkOptfJFbCQBQKcdeJvF2RjHSjSGYnRAxpnF8wbdmxUKjUuBd9hwmIqJxjmGYxr1wNDt9P1ZRi625jq4SKsXYVoaLaltx6GwDrpkXDyF8f8e54UQEaXDJVBM2HSmF2SrvNtVERES+hGGYxr1GOG9x1hxgwoEzdahvNUOtFGNqrfZRVhkA9PQxHg+umhOHhjYLdhVUe3sqREREbsMwTOPeNiyDGX0faJPUgWhc8gjsErDteFVXmcToVoYlScKHR8qwIDkScRGBckzZJyxPNyA8UI3NR8q8PRUiIiK3YRimcS9HTMFHWIMaaxDsElBi1yN/4RNIuPBWxIRrsfVYBdQKBSyjXBk+XtGM/KoWfGt2nMwz9y6NSoHLZkRja24l2s02b0+HiIjILRiGyS/kiCn4i/02aH/Xjg2qp/G7khkQQmDNVBN25FfDYrePurXa5qwyKBUCl00/f3sLD2b9zFi0mW3YdrzS21MhIiJyC4Zh8hsBAQEIDg7G2gRg+4lqnKxsxiVTo9FhsaO4rn1UZRKSJGHzkTIsS9NDFxLghll716JUHYyhAT010UREROMNwzD5lejoaMR0FkOrVuDl3WewKDUKoVpHPfFoWqsdOluP0oZ2bJg9fh6c602pELh8Zgy2n6hGU4fF29MhIiKSnUvbZAkhLgXwJAAlgBckSfp9v+PfAfBw18sWAHdLkpQl50SJ5GA0GnEiJwtXrVqE9w+V4qG1k7F6shGbjpSh0zLyMPxRVjk0KgXWTDW5Yba+4VuzYvHy7kJsPVaJa+fFe3s6srHa7Pgyrwrb8ipx6Gw9yhs7YLVJiAhSY0pMGJal6bFhdiyMYVpvT5WIiNxo2JVhIYQSwNMA1gGYCuAGIcTUfqedAXChJEkzAfwawHNyT5RIDjqdDseOHcN3L0hBp9WO1w+cxSXTHLW++87UjmgsSZKw9VgFVqQbEKpVu2O6PmF2QgRiw7X4/FiFt6ciC7tdwtsHi7HsD9tx138y8UVeJVINIdi4IBG3LUvGiokGVDR24IkteVjy+6/wyPvZqGzq8Pa0iYjITVxZGV4IoECSpNMAIIR4E8AGALndJ0iStKfX+fsAjJ/lIxpXjEYjDhw4gEnRoViWpsdre4vw+YMrAADSCEuGj5Y2oqyxA//vkklumKnvEELgkmnReOPAWbSZrQjSuPQXSj6puK4ND7x5GIfPNmBOYgR+feV0rJpshFIxcKOUU9UteGVPId44cBYfZ5fj8fXTcPXcuHGxqQoREZ3jSs1wHIDiXq9Lut4bzPcBfOrsgBDiDiFEhhAio7qajfzJ8/R6PSorK2GxWHDb0mRUNHVgR/6576I0gkT8+bEKKBUCF08xumOqPmXttGh0Wu345sT5++/t9uNVuOzvO1FQ1YK/fnsW3r/7AqyZanIahAFggiEEv9owHV/86EJMjg7Fj9/Jwi8/zBnT5ixEROR7XAnDzv5L4TQxCCFWwhGGH3Z2XJKk5yRJmi9J0nyDweD6LIlkolarER4ejtzcXKycZESyLggv7T7Tc/x4hfOtm535LKcCi1OjEBGkccdUfcqC5EhEBqnx2XlaKrHpcCl+8GoGEqOCsOWB5bhqjuvbZifrg/HmHUtw54Wp+M++s/jBqxnosLDvMhHReOFKGC4BkNDrdTyAAX2WhBAzAbwAYIMkSSMrviTyIJPJhH379kGhELj1gmQcPtsAVdfq4FfHq1wao6CqGaeqW7F22vjrLeyMSul4SPCrvCqYrefXyugHh0vw4FtHsDA5Cm/duQQJUUEjHkOpEHhk3RT89qoZ+PpENe7+TyY6rQzERETjgSth+CCAdCFEihBCA2AjgM29TxBCJAJ4H8DNkiSdlH+aRPIxGo04dOgQAOC6+QkIDVD19Bjelufa5hKfH3Ocd8lU/wjDgKNUornTij2narw9FZftOFmNn76TjSWpOrx82wKEBIyt3vnGRYn47VUzsP1ENX7yTvaIymqIiMg3DRuGJUmyArgPwOcA8gC8LUnSMSHEXUKIu7pOexSADsAzQogjQogMt82YaIy6O0oAQEiACtfNP/cXH4eLG1Db0jnsGJ8fq8DshAhEh/tP262laXoEa5Q9fxDwdScrm3H3fzKRbgrFs7fMg1atlGXcGxcl4uFLJ+OjrDI8vb1AljGJiMh7XNp0Q5KkLZIkTZQkaYIkSU90vfcvSZL+1fXrH0iSFClJ0uyuf+a7c9JEY2E0GnHq1Kme17csSer5tSQBXw/zkFhZQzuySxr9pkSim1atxEWTjfgitwL2UezW50mtnVbc899DCNSo8PJ3FyBM5tZ3d12YiqvmxOH/tp7El7nnxx8OiIjIOe5AR35Hp9OhpqYGnZ2OFeBkfXDPMX2IZti64e0nHMfXTB3/XST6WzPFhJoWM7JLG709lUFJkoSff3AUp6tb8PcbZrtl9V4Igd9dPQPTYsPw03ezUMU+xERE5y2GYfI7KpUKkZGROHz4cM97+pDujhACO05WwzJE+6ztx6uQEBWICYYQN8/U91w40QCFcP1BQ2/45Gg5Nh0pw4MXT8QFE/Ruu49WrcSTG+eg3WLDj9/J8vnVciIico5hmPxSdHQ09u7d2/O6u+ShpqUTzZ1WHCysc3pdh8WG3QW1WDnJ6JebL0QGazAnMRLbfTQM17Z04rEPj2FWfDjuuWiC2++XZgzBLy6fip35NXj9wFm334+IiOTHMEx+Sa/X91kZDuz3cNVXec7D3r7TtWi32LBysv+VSHRbNdmIo6WNPlka8PhHuWjqsOCP186CSumZ396+sygRF0zQ4Q+fHUd18/APXxIRkW9hGCa/pNfrkZvbs6M4lErHKm/3ZmSDlQF8faIaWrUCS1J1bp+jr1o5yfEHgeEeNPS0nfnV+CirDPetTMek6FCP3VcIgV9fOR2dFjt+80nu8BcQEZFPYRgmv2Q0GlFYWNjzWqUQUClET9A7XdOKMzWtfa6RJAlfHa/C0gl62dp0nY+mxIQiOkzrU3XDVpsdv/ooF4lRQbjrolSP33+CIQR3XTQBHx4pw56C86cPMxERMQyTn9LpdGhpaUF5eTkAQKlQwGqXcOOixJ5z+oe9U9WtOFvX5tclEoBjJXTlZAN2FdT4zG50/91/FvlVLfj55VMQoPLOH1TuuWgC4iMD8cSWPD5MR0R0HmEYJr+kUCgQFxeHrVu3AkDPdszL0vWIiwgEAHx1vG//2K+7Wqr5exgGHKUSLZ1WZAzyoKEnNbZZ8JcvTmJpmg6XTDV5bR5atRI/XTsJx8qasDlrwI71RETkoxiGyW/FxcXhm2++AQCoumqGJQn49gLHjnS7C2rR3GHpOf+r41WYZArtCcv+bGmaHhqlwidKJZ7feRqN7Rb8/LKpXu/wsX5mLKbHheFPn59Ah8Xm1bkQEZFrGIbJb8XFxSEjw7FzePfKsM0u4Zp58T3n7DjpqP9s7Wq3dtFkg+cn6oOCA1RYkBKJnfnerY+tazXj5d1ncPmMGEyNDfPqXABAoRD4n3VTUNrQjv/uZ6s1IqLzAcMw+a24uDgUFBTAZrNBqXD8q2C1S4iLCOzpFtG98nngTB0sNgkr0hmGuy1LM+BEZbNXW6w9u+MU2iw2PHhxutfm0N8FaXosSdXh2W9OcXWYiOg8wDBMfis0NBQBAQE4ePBgn5VhANi40FEq8d6hEkiShJ35NQhQKTAvKdJr8/U1y9Mdu7vt8lL3hOrmTry6pwgbZsUi3eS5VmquuH91GqqaO/FORrG3p0JERMNgGCa/lpycjM2bN0PZFYatXdswd+9IBwC55U3YVVCNhSlR0KqVqK2txY9//GPceeedKCgo8Mq8fcHUmDBEBWuwy0ulEi/tPoNOqw0PrPadVeFuS1J1mJcUiX99c9pnOm4QEZFzDMPk15KTk7Ft2zaoux6gs3atDGvVSqzp6kzw4q4zOFnZgmVpeuTn52PatGn47LPPkJmZiQULFuDo0aNem783KRQCF0zQYVdBDSTJs63EWjqt+M++Ilw6PRqphhCP3tsVQgjcvyoNpQ3t+OBwibenQ0REQ2AYJr+WmprqCLN2x+qdrVd/2PtWpgEA3j9UCgCYGxuENWvWYPLkybj++uuxfv16zJ07F9/+9rc9HgZ9xfJ0PaqaO3GyssWj933zwFk0d1hxx4oJHr3vSFw40YAZceF4dsdp9h0mIvJhDMPk18LCwhAaGoqczL0Azq0MA8DM+PCeX6uVAn97/CcIDAzERRdd1PP+0qVLUVlZiddff91jc/Yly7oeKNyZ77mtmS02O17adQYLU6IwOyHCY/cdKSEEfrA8BaerW/GNBz8fIiIaGYZh8nvp6enY9flmAIDNfq6+UwiBnyccxS7NAzihugG/M36En39rYp9etkqlEsuWLcMf/vAHj8/bF8RFBCLVEOzRFmtbjpajrLEDdyz3/LbLI7VuegyMoQF4eXeht6dCRESDYBgmv/fdecF4ZcIWnA64EYmvLgKy33YcyH4b36v7K+IVNVAIIDFc4PrAXZgu5fW5fvr06cjPz8fJkye9MHvvW56mx/4ztei0eqaN2Eu7C5FqCMaq82AnQI1KgVuWJGHHyWrkVzZ7ezpEROQEwzD5telSHm6JzEBiGKAQgKalFPjoAUcg/vRhKG19e+hqYMVq7Or7nkaDKVOm4KmnnvLk1H3GsnQDOix2ZBbVu/1eOaWNyCpuwM2Lk6BQeHe3OVfdsDARASoFXt5T6O2puEySJFQ3d+LQ2XocLKxDQVUzu2IQ0bil8vYEiLxpNXZBA2vfNy3tMH/8U6jNDXAWt8IxcIVv0qRJ+Pzzz90zSR+3ODUKSoXA3lO1uGCC3q33+s++IgSqlbh6bvzwJ/sIXUgArpwdh/cPleBn6yYjTKv29pQGdaq6Ba/tLcLWYxUoa+z3B0GlAksm6HDVnDhcPjMGaiXXUohofGAYJr/mLNgCgLqzAWKQhcdGDNzgITU1Fe+//z4qKythMpnknKLPC9WqMT0uHPtO17r1Pk0dFnx4pAzfmhWL8EDfDZTOfGdxIt7KKMaHR8pw8+Ikb09ngOrmTjzxSS4+zCqDWqnAinQDvr88FSn6IKiVCtS0dOJoSRO+zKvEg28dwf9tPYFfXjG1Tz9uIqLzFcMw+bVGhCLCWSAeJAhLALZh2YD3NRoNkpKS8MYbb+DBBx+UdY7ng8UpUXhp9xm0m20I1Cjdco/3M0vQbrHhJh8Mk8OZEReOqTFhePPAWZ8Lw5/lVODh97LRbrbhrgsn4AfLUqALCRhw3lVzgF9cPgXbT1ThT5+fwJ2vZeKKmTH4wzUzERzA/5QQ0fmLf89Ffm0blsHc78+EbRagrt35+W3QIkdMcXosJSXFj0sldLDYJBw+6566YUmS8N/9ZzErPhwzerW8O18IIXDDwgQcK2vC0ZJGb08HgOMz/esXJ3HXfzKRrAvClh8uw8OXTnYahLspFAKrp5jw0f3L8NO1k7DlaDmuemY3iuvaPDhzIiJ5MQyTX8sRU/AR1qABoZAANCAUb7ZdgA86Fg8IyWao8BlWDjpWYmIisrKy3Dxj3zQ/ORIKAbeVShw624D8qhbcuCjRLeN7woY5cdCqFXjj4FlvTwWSJOGxzcfw5LZ8XD03Dm/duQRpxoHlP4NRKxW4d2UaXvv+IlQ2deL6Z/fidLVnN14hIpILwzD5vRwxBU+K2/Er8f/wpLgdxRGLURp1wYCQ/BHWDLoqDAAxMTGora1FSYn/bb8bqlVjRlw49p2uc8v47x8qgVatwGUzYtwyvieEadW4fEYsNh8pQ2undfgL3ESSJPzvR7l4dW8R7liRij9fNwta9ehKW5am6fHG7Ythttrx7ef2cYWYiM5LDMNEg+gfkocKwoBjA474+Hh88sknHpqhb1mcqsOR4ga0m+XtN9xhseGjrDJcOi0aoT7cicEVNyxMQEunFVuOlnttDi/uOoN/7ynED5al4JF1k/tsIjMaU2PD8MYdi9FpseHWlw+goc0s00yJiDyDYZhIRgkJCfjqq6+8PQ2vWJyqg9lml71ueFteFZo6rLhm3vnTTm0w85IikaQLwqYjpV65/5e5lXhiSx7WTY/G/1w2ZcxBuNtEUyiev2U+Suracc9/D8HWa1tzIiJfxzBMJKP4+HjWDctcN/zeoRJEh2nd3sPYE4QQuHJ2HPacqkVFvz6+7lZS34YfvX0E02PD8ZfrZ8u+acmiVB1+c9V07DlViye/9M/dGIno/MQwTCSjmJgYFBUVwW73v926euqGz8hXN1zd3IlvTlbjyjlxUJ4nO84N58o5cZAk4KOsMo/d02qz44dvHgEk4JnvzHVb+7vr5yfgmrnx+Mf2AuzMr3bLPYiI5MYwTCSjkJAQaDQaZGZmensqXrE4VYcjZxvQYZGnbvjDI6Ww2SVcMzdOlvF8QYo+GLMSIvDBYc+VSjy1vQCZRfX4zVXTkRAV5NZ7/frKaUjVB+Ohd7PR3GFx672IiOTAMEwks7i4OHzzzTfenoZXdNcNHyqSp254c1YZZsSFI93ketuv88FVs2ORW96EExXOd0CUU35lM57eXoBvzYrFhtnu/0NFkEaF/7tuFiqbOvC7T4+7/X5ERGPFMEwkM5PJhAMHDnh7Gl4xLzkSQgAHCsdeKnG2tg3ZJY1YP+v8bac2mCtmxUKpEG5/kM5ul/A/HxxFcIAKj66f6tZ79TYnMRLfX5aC1/efxd5T7t2mm4horBiGiWQWGxuLnJwcb0/DK8K0akyODkOmDCvDn3S1HzufewsPRh8SgOXpenyUVQZJcl/nhbcyinGwsB4/v2wK9EPsLOcOP75kEhKiAvH45mOw2vyvhp6Izh8Mw0Qy8+eH6ABgflIkDhXVjzkAbTlajtkJEYiPdG+Nq7dcNiMGJfXtyCltcsv4TR0W/OnzE1iYEoVrvdCWTqtW4ueXTcWJyma8fsD7u+4REQ2GYZhIZqGhoRBCID8/39tT8Yr5yZFoNdtwfAz1sEW1rTha2ogrZo6/VeFul0w1QaUQ2JLjng04ntl+CvVtZjx6xVTZ+gmP1NppJlwwQYc/bz2J+lZuxkFEvolhmMgNTCYT9uzZ4+1peMW8pEgAGFOpRHeJxLpxWCLRLSJIgyUTdPj0aLnspRLFdW14afcZXDUnDtPjwmUdeySEEHhs/TQ0d1jwj68KvDYPIqKhMAwTuYHBYPDb9mpxEYGICdciYyxhOLsccxMjEBcRKOPMfM+66TEorG1DXrm8XSX+8sVJCAA/uWSSrOOOxqToUFwzNx7/2V/k8Y1GnGloM+OL3Eq8uOsMnvwyHy/uOoOvT1ShsZ1t4Ij8lcrbEyAaj4xGo98+RCeEwLykSGSMsqNEYU0rjpU14ZdXeK77gbdcMs2EX2w6ik9zyjE1NkyWMU9Vt+DDI6W4fXkqYn3kDxMPrE7HpiOleGp7Pn5z5QyP31+SJHx1vAr/3lOIXQU1cLYQr1IIrJpsxG1LU7Bkgs7jcyQi72EYJnIDo9Hot9syA46H6D7OLkdpQ/uIV3c/P1YBALh0erQ7puZT9CEBWJSiwydHy/H/1kyUpbb36a8KoFEpcPuKVBlmKI+EqCB8e0EC3jpYjDtXTHD7xh+95Vc249EPj2Hv6VrEhGtx/8o0LJ9oQLoxBCEBKjR1WHG8oglfn6jG+4dKsTV3Hy6aZMATV80Y938zQUQOLJMgcgODwYCKigpYLP75V6/zk6MAYFSrw1/kVmJabJjfBJHLZkTjdHUr8qtaxjxWYU0rNh0pxU2LkjzeSm04961MhxACz3ztmdphSZLwdkYx1j+1C3kVTfj1hmnY+dBK/L9LJmFBchQigjRQKRWICtbgggl6/M9lU7Dr4ZX4n8sm4+CZOlz61x349Kh7Hm4kIt/CMEzkBlqtFkFBQTh8+LC3p+IVk6NDEaxRIqNwZHXDtS2dyDxbjzVTTW6ame9ZO82xAv55TsWYx3pqewHUSgXuuNB3VoW7RYdrcd28eLx3qBRVze6tHZYkCb/dkoeH3s3G3MRIbP3RCty8JBkq5dD/ydOqlbhjxQR8+sMVSDOF4O7/HsIzXxe4tRc0EXkfwzCRm/jzTnQqpQJzEiNH/BDdtuNVkCTg4in+E4aNYVrMSojAl8erxjROeWM7Nh0uxQ0LE2EM1co0O3ndvjwVFpsd/95d6LZ72O0SfvbeUTy/8wxuWZKE176/aMSfR6IuCG/cvhjrZ8Xij5+dwJ+3nnTTbInIFzAME7mJwWDAoUOHvD0Nr5mXFIkTFU1o6nC9VOSL3ErEhmsxTaaHyc4XF082Iqu4YUwrpq/sKYJdkvD9ZSkyzkxeyfpgrJsejdf2FaGl0+qWe/x2Sx7eyijG/avS8L/fmgalYnR12Fq1Ek9+ezZuWJiAp7YX4B/b/LNvOJE/YBgmchODwYDc3FxvT8NrFiRHwS4Bh882uHR+h8WGnfnVuHiqyWubRHjL6q6V8O2jXB1uM1vxxoGzWDst2qMPp43GnSsmoLnDijfdsCvdCztP44VdZ/DdC5JleSBRoRB44soZuHpOHP78xUl8eKRUppkSkS9hGCZyE6PRiKKiIm9Pw2tmJ0ZACODwWddKJXbl16DDYverEoluU2JCERuuxZd5owvD7x0qRWO7xadXhbvNSojAopQovLy7EDa7fLW4uwtq8NstebhsRrSsu+4pFAK/v2YmFiRH4qF3s3G0pFGWcYnIdzAME7mJXq9HTU0NOjs7vT0VrwgJUGGiMdTlleEv8yoREqDC4lT/6/EqhMDqKaauPxDYRnSt3S7h5V1nMCs+vGf3P19329JklDa046sx1kl3q2jswA/fPIxUQwj+dO0sKEZZGjEYjUqBf940D/qQANz1n8wRlf4Qke9jGCZyE7VajfDwcGRkZHh7Kl4zOyECWSUNwz6Nb7dL+DKvChdOMkCj8s/fllZNMaLdYsPe07Ujuu6bk9U4XdOK7y1LOW/KSy6eYkJ0mBav7i0c81h2u4QfvXUEbWYb/nXTXAQHuKd9vj4kAP+4cQ4qmjrw6Cb/3FCHaLzyz//qEHmIP3eUABylEg1tFhTWtg153tHSRtS0dOLiKUYPzcz3LEnVIUijxLa8yhFd99/9Z6EPCcBlM2LcNDP5qZQK3LgoETvza3C6emz9lf974Cz2nq7FL6+YijRjqEwzdG5uYiTuX5WGTUfKsDmrzK336q+xzYKS+jZUNXew1RuRzLgDHZEb6fV6ZGdne3saXjM7IQIAcKS4Hin64EHP++ZkNYQAVqQbPDQz36NVK7EsTY9teVX49QbJpVXeisYOfHW8EndeOAHqYXro+pqNCxPw9235+O/+s6Peerukvg2/35KH5el6bFyQIPMMnbtvZRq+PlGNxzcfw4p0PSKCNG65j93u2EJ6c1YZ9p6uRXXzuXKrQLUS85MjsXZaNK6eG4cgDf9TTjQW59fvnkTnGYPBgLy8PG9Pw2smmkIRpFHiyDB1w9+crMaMuHDofGzXNE9bPcWI8sYOnKx0bbX0nYxi2CV4LAjKyRiqxaXTo/FORvGI66QBx8YaP//AUa7wu6tneKxERKVU4LdXzUBjuwV/+Oy4W+6x/XgVLn1yB37wagZ25ldjeZoeP79sCv547Uz877em4fr58ShtaMcvNuVgye++wvM7TsNis7tlLkT+gH+cJHIjg8GAXbt2eXsaXqNUCMyIC8eR4oZBz2lss+Dw2XrcuzLNcxPzUSsmOlbGd5ysxqToof/K32aX8ObBYixN0yFJN/iquy+7cVEiPs4ux+fHKrBhdtyIrt2WV4VvTlbj0SumIj7Ss+3kpsaG4XtLk/H8zjO4dl485iVFyTJuS6cVj27KwfuHS5FqCMaTG2fjshkxg676ZxbV4e/bCvDEljy8f7gUz3xn7pB/A0NEznFlmMiN9Ho9amtr/bajBADMSYxEbnnToKt/uwpqYJeAiyb5b4lEt5jwQKQbQ7Ajv3rYc3fmV6O0oR0bFyR6YGbusThFh/jIQLybWTKi6zqtNvzmk1ykGUNw85IkN81uaA9ePBGx4Vo8+uEx2GVoEVdc14ZrntmDD7PKcP+qNHz6w+XYMDtuyPKXeUlR+PdtC/Cvm+ahvLEd3/rHLnx1fGQ150TEMEzkViqVChEREcjMzPT2VLxmdkIELDYJx8qanB7/5mQVwrQqzIqP8OzEfNSKiQbsP1OHdvPQpQNvHihGVLAGl0w7f/syKxQC18yNx66CGpQ2tLt83b93F6Kwtg2/vGKq12qlgwNUeOjSyThW1oSPssf2MF1BVQuuemYPyhvb8e/bFuDHl0xCgErp0rVCCFw6PRqfPLAcyfpg3P5qJjcHIRohhmEiN/P3jhJzEiMAwGmphCRJ+OZkNZanG6A6zx4Ac5fl6XqYrXbsPzN4i7WGNjO2Ha/ElbPjXA5NvuraefGQJOB9F1eHG9sseGp7AVZNNuLCid7924RvzYrFtNgw/OnzE+i0jrzuGQDO1LTixuf3AZDw3t0XYPkoHyKNiwjEG3csxvykSDz41hF85OFuF0TnM/7Xh8jNDAYDsrKyvD0NrzGFaRETrnUahk9UNqOyqdProcaXLErRQaNSYMfJmkHP+eRoOSw2CVfPHVmdrS9KiArCklQd3j1U4lLLsBd2nUZzhxU/XTvJA7MbmkIh8Mi6KSipb8dre0e+22RtSydufnE/rHYJr9++GOmmsbWGCwlQ4ZXvLcSCpCj8+O0s7D01sp7VRP6KD9ARuZler/frjhKAo1TiSPHAbZm/OeGojV3BMNwjUKPEopQo7ByibnjT4VKkGUMwLTbMgzNzn2vnxePH72ThYGE9FqYM/jBaXasZL+06g8tnxGBKjG/87MvS9Viersc/vz6F7yxKQqDGtZV6s9WOu/97CNXNnXj7ziWYOMYg3E2rVuL5W+bj2n/twZ2vZeCTB5YjIWrsDxja7RJ2n6rBNyeqkV3aiJrmTgjh2IxkZnw4Vk4yYnGqTvbd/4g8gSvDRG5mMBhQWFjo7Wl41eyECBTXtaO2pe+DhN+crMbk6FBEh2u9NDPftCLdgPyqFpQ5qaMtrmvDwcJ6XDUn7rzZcW4462ZEI0ijxAeHh651fXbHKbRZbHjw4nQPzcw1P1ydjtpWM14/cNbla574JBcHztThj9fOxKyuftxyCQ9S46XvLgAA3P3fzFG1rutms0t4O6MYK/60HTe/eACv7iuCzS5hSkwYJkeHwWKz45U9Rbjxhf1Y+eev8V5mCWwyPFBI5EkMw0RuptPpUFtbi46ODm9PxWvObb7R0PNeu9mGjMJ6LE/Xe2dSPqx7pdzZ6nD3w1HfmhXr0Tm5U5BGhTVTTfg0p3zQfrm1LZ14dU8RNsyKHXM5gdzmJ0dhSaoOz+045VLw3H68Cq/sLcL3l6WMuKWcqxKigvDn62cjp7QJv/kkd1RjnKpuwZVP78ZD72ZDFxKAf9wwB9mPXYL37r4AT39nLp7+zly8f89SZD9+CZ7cOBuhWhV+/E4Wrv3XHhTWtMr8ExG5D8MwkZupVCpERkbi8OHD3p6K18yID4dSIXC41+YbmUX1MNvsuCCNYbi/iaYQRIdpB9QNS5KEDw6XYkFypCx/9e1L1s+MRUObBbvynddKv7q3CO0WG+5b5Zv9qO9flYbKps5h28TVtHTip+9mYXJ0KB661L11z2ummvCDZSn4z76z2HFy+HZ9vX2cXYb1/9iFkvo2PLlxNjbdcwHWz4qFVj2wDESrVmLD7DhsvncZ/nL9LJyqasFlf9+JL3PZ5o3ODwzDRB5gMpmwf/9+b0/Da4I0Kkw0hSKrpKHnvd2naqBSCCxMlmfDgvFECIFl6XrsPlXTp4ftsbImnKpuxZVzzv8H5/pbMdGAMK3KaReEdrMNr+4txMVTjEgz+taqcLclE3SYmxiBf359CtYhdoP7xQc5aOqw4m8bZ3ukE8hP1k7CBEMwHnn/KJo7LC5d8/LuM7jv9cOYGhOGLV39jl0pyVEoBK6eG4/Pf7QCacYQ3PFaBl7bN/IHC/tr6bQir7wJR0saUVzXxjIMkh0foCPyAH/vKAEAs+LD8dmxCkiSBCEE9pyqxeyECAQH8LchZy6YoMO7mSXIq2jCtNhwAMDmrDKoFAKXz4jx8uzkp1EpsG56DD45Wo4Oi63PCuS7mcWob7PgjhUTvDjDoQkhcPdFabj91Qx8fqwSl88c+P/RV8cr8dmxCvx07SRMjvbMA4BatRJ/um4Wrv3nHvzhs+P4zZUzhjz/hZ2n8ZtP8rB2mglPbpzjdCV4ODHhgXjzjsV44I3D+OWmHAgANy0e2eYopQ3teGP/WXx2rAIFVX23Jw9UK7EsXY+r5sRh7bRoKPnQHo0RV4aJPECv1yM3d3R1e+PFzPgINLRZcLauDU0dFhwtacAFE3TenpbPWtL12XS3x5IkCZ/mlGNpmh4RQRpvTs1t1s+KRUunFV+fqOp5z2aX8MKuM5idEIEFyZFenN3wVk02IjEqCC/tPjPgWIfFhsc2H0OaMQS3L0/16LzmJkbiliXJ+O/+s8gpbRz0vPcPleA3n+Th8hkxePrGuaMKwt2CNCo88515WDXZiF9synF5I5DGdgt+sekoLvzjdvzzm1Mwhgbgp2sn4ekb5+L5W+bj91fPwDXz4pBT2oh7/nsIq//8NXfdozFjGCbyAHaUAGYlOFY3jxQ3YP/pOtglYMkE1gsPJiY8EKn6YOzpCsPHyppQXNeOy2ZEe3lm7rM4NQr6EA029yqV+CK3EkW1bbhzRarPd89QKgS+e0EyMovqkdWvr/Yz2wtQXNeOX2+YDo3K8//p/dGaiYgM0uB/PzrmtJ/zvtO1eOjdbFwwQYe/fHuWLJvgaFQKPPOduViYEoWfvpuN7F5lUs7syq/Bmr98g9f3n8UNCxOx46GVeP32xbh3ZRounxmDNVNN2LgwEb+5cgZ2PbwK/7ppLlRKBb737ww8+OZhtJmtI5qf1WZHcV0bMovqkV3SgKpm/33I2d/x7yeJPECn06Gurg7t7e0IDAz09nS8YqIpFAEqBbJLGmGXJASoFJibFOHtafm0JRN0+PBIGaw2Oz7NKYdSIbBm6vgNwyqlo1TincxitJmtCNKo8J99RYgN1+KSaefHz33d/Hj85YuTeHn3Gfxt4xwAQEVjB57dcRrfmhXbs+LvaeGBajy0dhJ+9v5RbM4q69PFoqqpA/e9fhiJuiA8e/M8WWuZtWol/vmdufjWU7txx6uZ2Hz/UhhD+7ZSlCQJL+8uxG8+yUWaMQQv3roAM+LDhxxXqRC4dHoMVk024ZmvC/Dktnwcr2jG87fMH/LhUrtdwjf51Xj7YDF25degubNvgDaFBeCSqdH49oIETI8beg40fnBlmMgDlEoloqKikJmZ6e2peI1aqcD0uHBklzRg76laLEiOOu+3Ena3Cybo0dJpRXZpIz49WoHFqVGICh6fJRLd1k2PRofFjh0nq3G6ugW7Cmpww8LE86YuNFSrxnXz4/HJ0XJUNjlWGv/25UnYJcnru+ZdNz8B0+PC8MfPzm0fbbNLuO/1w2jttOJfN81DqFYt+311IQF4/pb5aGg34yfvZPdZmZYkCX/94iR+9XEu1kw14YN7lg4bhHvTqBR48OKJeOW2hShraMe3n92LM4O0dcsorMPl/9iF214+iIOF9bh8Zgx+f/UMvHzbArxwy3w8esVUzEuKxNsZxbjiH7twx6sZKK5rG3YO7WYbDp2tx9ZjFfgytxJHihvG1NuZPM+llWEhxKUAngSgBPCCJEm/73dcdB2/DEAbgO9KknRI5rkSnde6O0osW7bM21Pxmpnx4fjv/rMwW+346drx0yfXXRanOjptvLKnEKdrWvG9ZSlenpH7LUyJQmSQGp/lVOBgYT1UCoFvL0zw9rRG5NYlyXh5dyHePliMdTOi8XZGMW69INnr7fCUCoGH1k7GLS8dwJsHHHN6fudpHCisw1+unyXbLnjOTI0Nw88vm4JffngMr+4twq0XJAMA/vZlPv7+VQGunx+P3189c9Q72K2YaMAbdyzGzS8ewLef3Yv37r6g5/O22Oz43ZbjeGn3GcSGa/Hn62Zh/axYp+Uq31uWgsZ2C17ZU4hnvzmFtX/bgcfXT8P1C/p+B602O7bkVOCtg2dx4EwdLLa+pScapQJLJuiwcUECLhnmIb/WTitaOq0ICVDxgWIvGfZTF0IoATwNYA2AEgAHhRCbJUnq/TTQOgDpXf8sAvDPrv8loi7sKAFcJu3C9xV/RGxADawZcYDucWDm9d6els/ShQTgrshM3JT7AP4aUANpdzwQ9Ni4/sxUSgV+GpuNC/P+iRjU4p4gI3RnfnNe/czJ+mD8NCYL1+z8IUw7a7BLo0Oo6dcApnl7aliersePTEdwydYHIG2txXpJB3XiHbhqzmVuv/dNi5Ow7XgVsj99HpY9H0DVUobr7DoYJ9yNG66+bMxbOU+LDcebdyzGtf/cg9ee+xN+pnkboqkUDQoDajquxa1LvoOH101GkGbo6BMeqMYDq9Nxzbx4PPRuFh56LxvKnHdwdcOLEI2l6AiKwV/sG/Fcw3wkRAXie8tSMC8xEtHhWkgSUN7YgcyiOnycXY67/3sIk6ND8fdpBZiY8xegsQRSeByOpP8Qz9bPw4HCOtS1mnvuHR2mxYKUKFw5OxYXdX4N5fZfA40lQHg8zBf9Ann6S1Hd3Am7JEEfGoAp0WEIPP4esO1XPedh9aN9/33Jfnvw486OAX3fS78EyN/q/HVg10Ot7XUABIARtL17fPCHOT1NOCuk73OCEEsAPC5J0tqu148AgCRJv+t1zrMAvpYk6Y2u1ycAXCRJUvlg486fP1/KyMgY+08wQn/605/Q1jb8X3sQyS03NxenT5/GwYMHvT0V78h+G/bND0Bh7bXFsDoQWP/38yroeFT22zB/cB80Uq9trMf7Z5b9Nmwf3g+lrdfDTOfbz5z9Nqwf3g+VL/4M2W/D9uEDUNrO/XsoqQIhvuWZuTUe+C80Wx5EIM4FQEkdCCHjZ1Pw5UuI3fkwgsS5e1iVWqg2/GPE97Da7PjoP09i7enf9hmvAwE4sfA3mHHp7YOGeJtdwidHy3H442fxU/Mzfa5vkzT4veoetE66GhNNIQjRqtDYbsGJimbsyq/B0vav8EfNi9Cis881P7P8AJvt5/528SrVbvxe/QICev0e0SkC8FTw/fhCdSEuxy7c2fi3Pr+HWBVabJ3wP1CrFFh58jd9vqd2oQKEgMJ+rie1BEfMHez1mHgwEAshMiVJmu/0mAth+FoAl0qS9IOu1zcDWCRJ0n29zvkYwO8lSdrV9XobgIclSRo07TIMk7+prq7GW2+9herqke0ENW78dTrQWDzw/fAE4Ec5np/P+cAfP7Px8DP78s/g7bl54v5y32OM49n/Mg2KpoE7E0rhCRBOrrfY7LD83zQEtTvZgCYoFidv2AuFEChvbMeiDy9EuLliwHk1SiMeSXoDTxTeAKO9asDxErujk0+8wvmOjx5zHoXh6wCs7ReGF0qSdH+vcz4B8Lt+YfghSZIy+411B4A7ul5OAnBidD/SmOkBePkbMC7wc5SHX3yO82IU8wY7lllul+PJwnH3OXrgM3PGq5+jl35mWfnyz+DtuY3y/iP6Tsr9M451vNFc7+o1w5031HFf4OF/H5IkSTI4O+BKpXYJgN6V4/EA+v9xxZVzIEnScwCec+GebiWEyBjsTwfkOn6O8uDnKA9+jvLg5ygPfo7y4WcpD36Og3OltdpBAOlCiBQhhAbARgCb+52zGcAtwmExgMah6oWJiIiIiHzBsCvDkiRZhRD3AfgcjtZqL0mSdEwIcVfX8X8B2AJHW7UCOFqr3ea+KRMRERERycOlhnaSJG2BI/D2fu9fvX4tAbhX3qm5lddLNcYJfo7y4OcoD36O8uDnKA9+jvLhZykPfo6DGPYBOiIiIiKi8YrbMRMRERGR3xq3YVgIcakQ4oQQokAI8TMnx4UQ4u9dx7OFEHO9MU9f58LneJEQolEIcaTrn0e9MU9fJ4R4SQhRJYRw2pSS30fXuPA58vvoAiFEghBiuxAiTwhxTAjxQyfn8Ds5DBc/R34nhyGE0AohDgghsro+x/91cg6/jy5w8bPkd7KfcbkJNreQloeLnyMA7JQk6QqPT/D88m8ATwF4dZDj/D665t8Y+nME+H10hRXAjyVJOiSECAWQKYT4gr9HjpgrnyPA7+RwOgGskiSpRQihBrBLCPGpJEn7ep3D76NrXPksAX4n+xivK8MLARRIknRakiQzgDcBbOh3zgYAr0oO+wBECCFiPD1RH+fK50gukCRpB4C6IU7h99EFLnyO5AJJksolSTrU9etmAHkA4vqdxu/kMFz8HGkYXd+xlq6X6q5/+j/QxO+jC1z8LKmf8RqG4wD03juxBAN/g3LlHH/n6me0pOuvZD4VQkzzzNTGHX4f5cPv4wgIIZIBzAGwv98hfidHYIjPEeB3clhCCKUQ4giAKgBfSJLE7+MoufBZAvxO9jFew7Bw8l7/Pxm5co6/c+UzOgTHFoezAPwDwCZ3T2qc4vdRHvw+joAQIgTAewAelCSpqf9hJ5fwO+nEMJ8jv5MukCTJJknSbDh2sF0ohJje7xR+H13kwmfJ72Q/4zUMy7aFtJ8b9jOSJKmp+69kuvpRq4UQes9Ncdzg91EG/D66rque8D0A/5Uk6X0np/A76YLhPkd+J0dGkqQGAF8DuLTfIX4fR2iwz5LfyYHGaxjmFtLyGPZzFEJECyFE168XwvGdqvX4TM9//D7KgN9H13R9Ri8CyJMk6S+DnMbv5DBc+Rz5nRyeEMIghIjo+nUggIsBHO93Gr+PLnDls+R3cqBx2U2CW0jLw8XP8VoAdwshrADaAWyUuJPLAEKINwBcBEAvhCgB8BgcDzbw+zgCLnyO/D66ZimAmwEc7aotBID/AZAI8Ds5Aq58jvxODi8GwCtdHYwUAN6WJOlj/jd7VFz5LPmd7Ic70BERERGR3xqvZRJERERERMNiGCYiIiIiv8UwTERERER+i2GYiIiIiPwWwzARERER+S2GYSIiIiLyWwzDREREROS3GIaJiIiIyG/9f13v4Yy2E5nmAAAAAElFTkSuQmCC\n", 251 "text/plain": [ 252 "<Figure size 864x648 with 1 Axes>" 253 ] 254 }, 255 "metadata": { 256 "needs_background": "light" 257 }, 258 "output_type": "display_data" 259 } 260 ], 261 "source": [ 262 "plt.rcParams[\"figure.figsize\"] = (12,9)\n", 263 "fig = plt.figure()\n", 264 "plt.plot(res[:,0],res[:,1])\n", 265 "x_grid = np.linspace(0, 1., 1000)\n", 266 "plt.plot(x_grid, (1 - x_grid + 0.05 * np.cos(11 * np.pi * x_grid)), 'k-', linewidth=1)\n", 267 "plt.plot(b_points[:, 0], b_points[:, 1], 'o')\n", 268 "plt.fill_between(x_grid, -1., (1 - x_grid + 0.05 * np.cos(11 * np.pi * x_grid)), color='gray')\n", 269 "plt.ylim(0, None);" 270 ] 271 }, 272 { 273 "cell_type": "markdown", 274 "id": "d39bab91", 275 "metadata": {}, 276 "source": [ 277 "We can clearly see the dissipative behaviour of the inelastic collisions at play. Indeed, despite the fact that we asked to integrate the system up to $t=10$, the integration was stopped at" 278 ] 279 }, 280 { 281 "cell_type": "code", 282 "execution_count": 8, 283 "id": "f75d9d74", 284 "metadata": {}, 285 "outputs": [ 286 { 287 "data": { 288 "text/plain": [ 289 "4.739019126545039" 290 ] 291 }, 292 "execution_count": 8, 293 "metadata": {}, 294 "output_type": "execute_result" 295 } 296 ], 297 "source": [ 298 "ta.time" 299 ] 300 }, 301 { 302 "cell_type": "markdown", 303 "id": "6a619f29", 304 "metadata": {}, 305 "source": [ 306 "instead, and the outcome ``oc`` informs us that the integration was stopped because of a stopping terminal event:" 307 ] 308 }, 309 { 310 "cell_type": "code", 311 "execution_count": 9, 312 "id": "935109c2", 313 "metadata": {}, 314 "outputs": [ 315 { 316 "data": { 317 "text/plain": [ 318 "<taylor_outcome.???: -2>" 319 ] 320 }, 321 "execution_count": 9, 322 "metadata": {}, 323 "output_type": "execute_result" 324 } 325 ], 326 "source": [ 327 "oc" 328 ] 329 }, 330 { 331 "cell_type": "markdown", 332 "id": "09145c42", 333 "metadata": {}, 334 "source": [ 335 "A value of ``oc`` negative and greater than ``taylor_outcome.success`` indicates a stopping terminal event. The index of the stopping terminal event is ``-oc - 1``:" 336 ] 337 }, 338 { 339 "cell_type": "code", 340 "execution_count": 10, 341 "id": "ba9dc815", 342 "metadata": {}, 343 "outputs": [ 344 { 345 "name": "stdout", 346 "output_type": "stream", 347 "text": [ 348 "Stopped by terminal event? True\n", 349 "Index of the terminal event: 1\n" 350 ] 351 } 352 ], 353 "source": [ 354 "print(\"Stopped by terminal event? {}\".format(oc > hy.taylor_outcome.success and int(oc) < 0))\n", 355 "print(\"Index of the terminal event: {}\".format(-int(oc)-1))" 356 ] 357 } 358 ], 359 "metadata": { 360 "kernelspec": { 361 "display_name": "Python 3 (ipykernel)", 362 "language": "python", 363 "name": "python3" 364 }, 365 "language_info": { 366 "codemirror_mode": { 367 "name": "ipython", 368 "version": 3 369 }, 370 "file_extension": ".py", 371 "mimetype": "text/x-python", 372 "name": "python", 373 "nbconvert_exporter": "python", 374 "pygments_lexer": "ipython3", 375 "version": "3.8.10" 376 } 377 }, 378 "nbformat": 4, 379 "nbformat_minor": 5 380} 381